PINNForge: Execution-Grounded Evolutionary Design of Physics-Informed Neural Networks for PDE Solving via Large Language Models

# PINNForge: Execution-Grounded Evolutionary Design of Physics-Informed Neural Networks for PDE Solving via Large Language Models
Source: [https://arxiv.org/html/2609.23023](https://arxiv.org/html/2609.23023)
Xu YangAffiliation:East China Normal UniversityJun ZhangAffiliation:Hanyang UniversityXiaolong WangAffiliation:East China Normal UniversityJing XuAffiliation:Nankai UniversityKeqian LiAffiliation:East China Normal University

###### Abstract

Physics\-informed neural networks \(PINNs\) require coordinated choices over network representation, sampling, loss construction, and optimization, while effective configurations often vary substantially across partial differential equations \(PDEs\)\. Existing automated PINN design methods can search candidate configurations, but information revealed during actual training is still used mainly for evaluation rather than to improve subsequent design, leading to repeated trial\-and\-error and inefficient use of training budget\. We propose PINNsForge, an LLM\-driven evolutionary framework for execution\-feedback\-based automated PINN design\. PINNsForge generates diverse candidate configurations from PDE\-related prior knowledge, evaluates them through actual training, and feeds high\-performing designs together with accumulated execution evidence back to the LLM\. Guided by observed optimization behavior, the LLM then refines, recombines, and explores coupled PINN design components, forming a continual cycle of generation, execution, feedback, and evolution\. Unlike one\-shot search or evaluation\-only feedback, PINNsForge progressively converts training experience into improved design decisions for the target PDE\. Across 25 PDE benchmarks, PINNsForge achieves the lowest mean MSE on 24 tasks compared with RoPINN, PINNsFormer, and PINNsAgent\. Ablation studies further confirm the importance of the PDE knowledge base, execution feedback, and evolutionary search: removing these components increases the mean MSE to 3\.74×\\times, 12\.10×\\times, and 10\.10×\\timesthat of the full PINNsForge, respectively\.

## 1Introduction

Physics\-informed neural networks \(PINNs\) provide a unified neural numerical paradigm for complex physical systems by directly embedding partial differential equations \(PDEs\), together with their initial and boundary conditions, into neural network training\([Raissi et al\., 2019](https://arxiv.org/html/2609.23023#bib.bib27);[Toscano et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib34);[Bai et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib2);[Michaloglou et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib24)\)\. However,reliably solving a PDE with a PINN is often substantially more difficult than constructing the PINN itself\([Krishnapriyan et al\., 2021](https://arxiv.org/html/2609.23023#bib.bib16);[Rathore et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib28)\)\. PINN performance is jointly affected by input representation, network architecture, collocation\-point sampling, physical\-constraint handling, loss construction and weighting, optimizers, and the training procedure, with strong interactions among these factors\([Toscano et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib34);[Michaloglou et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib24);[Gao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib8);[Zhang and Ding, 2025](https://arxiv.org/html/2609.23023#bib.bib59);[Cao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib4);[Zhang and Deng, 2025](https://arxiv.org/html/2609.23023#bib.bib60)\)\. Prior studies further show that PINNs can suffer from gradient imbalance\([Wang et al\., 2021a](https://arxiv.org/html/2609.23023#bib.bib39)\), inconsistent convergence rates across loss components\([Wang et al\., 2022](https://arxiv.org/html/2609.23023#bib.bib41)\), insufficient local sampling\([Luo et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib22)\), and long\-horizon error propagation\([Lin and Chen, 2025](https://arxiv.org/html/2609.23023#bib.bib17)\)\.

More importantly,the dominant training bottlenecks differ across PDEs\.Multiscale problems may be limited primarily by the ability to represent high\-frequency components, complex geometries may depend more strongly on boundary sampling and constraint handling, long\-time dynamical systems may be limited by temporal causality and optimization stability, while high\-dimensional PDEs and inverse problems may exhibit entirely different training behavior\([Toscano et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib34);[Bai et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib2);[Thawon et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib33);[Michaloglou et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib24);[Wang and Hu, 2025](https://arxiv.org/html/2609.23023#bib.bib37);[Francis et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib6);[Francis et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib7);[Yang et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib50);[Liu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib21);[Ren et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib29)\)\. Consequently, a network architecture or training technique that is effective for one class of PDEs cannot be assumed to transfer directly to another\. An unsuitable configuration may not only reduce accuracy, but can also lead to persistent PDE\-residual stagnation, inadequate satisfaction of initial or boundary conditions, severe imbalance among loss terms, or even failure of the overall solution process to converge\. The central challenge across heterogeneous PDEs is therefore not the absence of a single stronger PINN method, but rather how to identify an effective combination of coupled design choices spanning network representation, sampling, physical constraints, loss construction, and optimization strategies\.

A large body of work has proposed effective techniques for specific training difficulties, including adaptive loss balancing, residual\-driven sampling, causal training, and new network representations\([Gao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib8);[Cao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib4);[Luo et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib22);[Wang and Hu, 2025](https://arxiv.org/html/2609.23023#bib.bib37);[Lin and Chen, 2025](https://arxiv.org/html/2609.23023#bib.bib17);[Zhang and Ding, 2025](https://arxiv.org/html/2609.23023#bib.bib59);[Zhang and Deng, 2025](https://arxiv.org/html/2609.23023#bib.bib60);[Yang et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib50)\)\. However, these methods typically assume in advance which problem should be addressed\. For a new PDE, it is generally unknown whether representation capacity, sampling, physical constraints, or optimization will become the dominant bottleneck, and it is even harder to determine beforehand whether multiple interdependent components must be adjusted jointly\. Thus, continually adding new PINN techniques does not resolve a more fundamental question:who determines which mechanisms are needed for the current PDE, and how should these mechanisms be combined?

This challenge has motivated the development of automated PINN design\([Zhou et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib63);[Liu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib21);[Ren et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib29)\)\. Auto\-PINN uses neural architecture search and hyperparameter optimization to automatically select configurations such as network depth, width, and activation functions\([Wang et al\., 2023](https://arxiv.org/html/2609.23023#bib.bib43)\)\. PINNsAgent further combines LLMs, Physics\-Guided Knowledge Replay, and Memory Tree Reasoning to search for PINN architectures and training configurations across different PDEs\([Wuwu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib49)\)\. Work such as Lang\-PINN also shows that LLMs can participate in PDE understanding, code generation, and diagnosis of execution results\([He et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib10)\)\. These methods substantially reduce the reliance of PINN design on manual trial and error and indicate that LLMs are increasingly capable of handling complex PINN design tasks\.

However, we argue that a key gap remains in current automated PINN design:the rich information produced during real PINN training has not yet been fully exploited to guide subsequent algorithm design and search\.A candidate training run provides far more than a final performance score\. Signals such as the PDE residual, initial\- and boundary\-condition errors, the evolution of different loss components, convergence stagnation, and numerical stability can reveal why a configuration works, where its limitations arise, and what specific design requirements the current PDE imposes\. Existing automated methods, however, tend to focus more on candidate generation and performance ranking, with limited use of the specific problems and their causes exposed during training\. The problem is therefore not merely how to automatically generate and evaluate more PINN configurations, buthow to extract information from real training that can influence the next design round, allowing the search to actively change direction according to the difficulties actually exhibited by the current PDE\.

This challenge is further amplified by the strong coupling among components in the PINN design space\. Performance bottlenecks observed during real training often cannot be resolved by modifying a single hyperparameter; they may jointly involve network architecture, activation functions, loss design, sampling strategy, and optimization\. Consequently, independent parameter adjustments or predefined local search operations are often insufficient to capture interactions among components or to produce coordinated modifications in response to the observed training behavior\.An automated design system therefore needs more than the ability to search additional candidates: it must understand the relationships among candidate configurations, training behavior, and failure causes, and use them to decide which mechanisms should be inherited, which components should be adjusted jointly, which existing designs should be combined, and when the search should shift toward new directions\.

This observation naturally connects to recent work on LLM\-driven automated algorithm discovery\. Studies such as Evolution of Heuristics \(EoH\), FunSearch, ReEvo, and LLaMEA show that LLMs can be embedded in a “generate—execute—evaluate—evolve” loop, where algorithms are iteratively modified and recombined at a semantic level based on existing candidates and their real execution results\([Liu et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib18);[Romera\-Paredes et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib30);[Ye et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib52);[van Stein and Bäck, 2025](https://arxiv.org/html/2609.23023#bib.bib35);[Wu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib48);[Yao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib51);[Zhong et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib62);[Yin et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib53);[Yu and Liu, 2025](https://arxiv.org/html/2609.23023#bib.bib55);[Liu et al\., 2026b](https://arxiv.org/html/2609.23023#bib.bib19);[Liu et al\., 2026a](https://arxiv.org/html/2609.23023#bib.bib20);[Yu and Liu, 2026](https://arxiv.org/html/2609.23023#bib.bib56);[Qin et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib26);[Wang et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib38)\)\. In contrast, the limited prior work applying evolutionary computation to automated PINN design mainly relies on predefined selection, crossover, and mutation operators, using final performance as fitness to search network architectures or hyperparameters\([Kaplarević\-Mališić et al\., 2023](https://arxiv.org/html/2609.23023#bib.bib13);[Zhang and Yang, 2024](https://arxiv.org/html/2609.23023#bib.bib58);[Carrillo et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib5)\)\. This suggests a largely unexplored direction: can EoH\-style LLM\-driven evolutionary search be introduced into PINN design so that training behavior exposed during real PDE solving directly participates in designing the next generation of solvers?

Building on this idea, we proposePINNsForge, anexecution\-feedback\-driven evolutionary search frameworkover a unified PINN design space\. PINNsForge turns each PINN training run from a mere candidate evaluation into part of the subsequent algorithm\-design process\. The system first generates candidates with different design preferences based on PDE characteristics and relevant PINN knowledge\. After real training, it extracts execution information including performance, PDE residuals, physical\-constraint satisfaction, and optimization behavior, and feeds this information back to the LLM together with high\-performing candidates from the search history\. The LLM then performs cross\-component modification and recombination over representation, network architecture, sampling, constraints, loss design, and optimization through complementary strategies including conservative refinement, physics\- and constraint\-guided correction, architecture and optimization adjustment, multi\-parent synthesis, and novelty exploration\.

PINNsForge therefore does not aim to propose another fixed PINN architecture or training technique, but instead addresses a higher\-level automated algorithm\-design question:

Can an automated PINN designer learn from the actual performance of previous solvers and continually evolve solvers that are better suited to the target PDE?

PINNsForge unifies algorithm design, real PDE solving, training\-behavior analysis, and subsequent search within a continuously iterated closed loop\. Real training results are used not only to answer “which candidate is better?” but also to determine “how should the next round change?” As a result, the LLM is no longer merely a one\-shot configuration generator, but becomes a search module continually driven by real PDE\-solution outcomes, extending automated PINN design from staticgenerate\-and\-evaluatetoward a closed\-loopdesign\-execute\-diagnose\-evolveprocess for algorithm discovery\.

The main contributions of this work are as follows:

1. 1\.We introduce an automated PINN algorithm\-design framework driven by real execution feedback\.PINNsForge extends automated PINN design from one\-shot architecture or hyperparameter selection to multi\-generation closed\-loop search, jointly considering input representation, network architecture, sampling, physical constraints, loss construction, optimization strategy, and training control within a unified structured design space\. Unlike methods that generate candidates primarily from static PDE descriptions, PINNsForge grounds subsequent search continually in the candidates’ actual PDE\-solving behavior\.
2. 2\.We introduce EoH\-style LLM\-driven evolutionary search into automated PINN design\.PINNsForge is the first work to apply the Evolution of Heuristics paradigm to complete PINN design\. Unlike conventional evolutionary search based on predefined mutation and crossover operations, PINNsForge jointly provides the LLM with PDE priors, high\-performing historical candidates, and performance, physical\-constraint, and optimization information obtained from real training\. This enables conditional inheritance, modification, and recombination of multiple coupled PINN components according to actual candidate behavior, while allowing the search to explore new solution strategies when existing designs are insufficient\.
3. 3\.We systematically evaluate execution\-feedback\-driven evolutionary search on a large and heterogeneous PDE suite\.We evaluate on 25 PDEs spanning nonlinear, multiscale, high\-dimensional, complex\-geometry, long\-time dynamical, fluid\-system, and inverse problems\. PINNsForge achieves the lowest mean MSE on 24 tasks in a unified comparison with RoPINN, PINNsFormer, and PINNsAgent\. In addition, component ablations on 12 representative PDEs, a budget\-matched non\-evolutionary search control, multi\-generation search trajectories, and low\-/high\-fidelity ranking analyses provide complementary evidence for the roles of PDE knowledge, real execution feedback, parent\-driven joint evolution, and multi\-fidelity validation in final performance\.

## 2Related Work

### 2\.1PINN Optimization and Failure Modes

Physics\-informed neural networks incorporate the PDE residual together with physical constraints such as initial and boundary conditions into the training objective, enabling neural networks to approximate PDE solutions without requiring large amounts of supervised data\. PINN training, however, is not a simple supervised\-learning problem\. Different physical constraints and loss components can have substantially different scales and optimization dynamics, leading to gradient imbalance, inconsistent convergence rates, and training stiffness\.

[Wang et al\. \(2021a\)](https://arxiv.org/html/2609.23023#bib.bib39)analyzed pathological behavior in PINN training from the perspective of gradient flow, showing that gradient imbalance across loss components can substantially affect optimization and motivating an adaptive weighting strategy\.[Wang et al\. \(2022\)](https://arxiv.org/html/2609.23023#bib.bib41)further analyzed the differing convergence rates of PINN loss components through the Neural Tangent Kernel \(NTK\), providing theoretical insight into understanding and balancing multi\-objective physics losses\.

Beyond loss balancing, the distribution of collocation points directly determines the physical supervision received by a PINN in different regions\.[Wu et al\. \(2023\)](https://arxiv.org/html/2609.23023#bib.bib46)systematically studied non\-adaptive and residual\-driven adaptive sampling and proposed methods such as RAD and RAR\-D, which allocate more collocation points to difficult regions according to the current PDE residual\. For time\-dependent PDEs,[Wang et al\. \(2024\)](https://arxiv.org/html/2609.23023#bib.bib42)further showed that ignoring temporal causality in the physical system can cause errors to propagate and accumulate over long horizons, and used causal training to alleviate optimization difficulties in complex dynamical systems\.

Beyond gradients and sampling, prior work has characterized PINN training difficulties from broader perspectives\.[Krishnapriyan et al\. \(2021\)](https://arxiv.org/html/2609.23023#bib.bib16)systematically analyzed failure modes in convection, reaction, and diffusion problems and showed that PDE regularization terms can induce severe optimization pathologies\.[Rathore et al\. \(2024\)](https://arxiv.org/html/2609.23023#bib.bib28)further studied ill\-conditioning induced by differential operators and the effect of optimizer choice from a loss\-landscape perspective\.

At the level of individual design components,[Jagtap et al\. \(2020\)](https://arxiv.org/html/2609.23023#bib.bib11)studied adaptive activation functions;[Wang et al\. \(2021b\)](https://arxiv.org/html/2609.23023#bib.bib40)improved multiscale representation using Fourier features and spectral\-bias considerations;[McClenny and Braga\-Neto \(2023\)](https://arxiv.org/html/2609.23023#bib.bib23)introduced adaptive loss balancing through trainable weights;[Sukumar and Srivastava \(2022\)](https://arxiv.org/html/2609.23023#bib.bib31)studied exact imposition of boundary conditions;[Yu et al\. \(2022\)](https://arxiv.org/html/2609.23023#bib.bib54)incorporated gradient information from the PDE residual into training; and XPINNs and hp\-VPINNs extended PINN designs through space–time domain decomposition and variational formulations\([Jagtap and Karniadakis, 2020](https://arxiv.org/html/2609.23023#bib.bib12);[Kharazmi et al\., 2021](https://arxiv.org/html/2609.23023#bib.bib14)\)\.

Together, these studies show that PINN performance can be jointly affected by network representation, sampling, physical constraints, loss balancing, and optimization\. However, these improvements are typically designed for specific training difficulties, whereas the dominant failure modes differ across PDEs\. Multiscale problems may be limited primarily by representation capacity, complex\-boundary problems may depend more strongly on constraint handling and sampling, and long\-time dynamical problems may be affected by temporal causality and optimization stability\. This motivates a natural question: rather than fixing a particular improvement strategy in advance, can a more suitable combination of designs be automatically discovered from the characteristics of the target PDE and its actual training behavior?

### 2\.2Automated PINN Design

As the number of tunable PINN design factors has increased, automated network design and hyperparameter search have been used to reduce the cost of manual experimentation\. Auto\-PINN introduces neural architecture search and automated hyperparameter optimization into PINNs, automatically exploring architectural choices such as network depth, width, and activation functions\([Wang et al\., 2023](https://arxiv.org/html/2609.23023#bib.bib43)\), showing that PINN architectures can be determined through systematic search rather than relying entirely on manual experience\. Meanwhile, PINNacle provides a large\-scale benchmark covering diverse PDE types and typical PINN failure scenarios, establishing a systematic foundation for unified comparison and cross\-problem evaluation of PINN methods\([Hao et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib9)\)\.

Beyond Auto\-PINN, prior work has explored automated PINN design through neural architecture search and evolutionary computation\. NAS\-PINN uses differentiable neural architecture search to automatically select the number of hidden layers, neurons per layer, and residual connections\([Wang and Zhong, 2024](https://arxiv.org/html/2609.23023#bib.bib44)\)\.[Kaplarević\-Mališić et al\. \(2023\)](https://arxiv.org/html/2609.23023#bib.bib13)proposed PINN/GA, which jointly searches network architectures and hyperparameters using a genetic algorithm\.[Zhang and Yang \(2024\)](https://arxiv.org/html/2609.23023#bib.bib58)further used evolutionary computation to jointly explore network topology, shortcut connections, and parameterized activation functions, while MOPINNs applies evolutionary multi\-objective optimization to loss weights and network hyperparameters\([Carrillo et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib5)\)\.

These studies demonstrate the feasibility of automated search and evolutionary computation for PINN design, but their optimization targets are concentrated primarily on network topology, activation functions, hyperparameters, or loss trade\-offs, and their search operations mainly rely on predefined encodings and evolutionary operators\. In contrast, PINNsForge focuses on using performance, physical, and optimization diagnostics produced by real training as context, allowing an LLM to semantically inherit, modify, and recombine multiple design components over multiple search generations\.

Recent studies have also begun to use large language models \(LLMs\) for automated PINN design\. PINNsAgent combines Physics\-Guided Knowledge Replay with Memory Tree Reasoning to retrieve prior PDE\-solving experience and search candidate configurations, enabling automated exploration of PINN architectures and hyperparameters across different PDEs\([Wuwu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib49)\)\. Compared with conventional numerical hyperparameter optimization, such approaches allow LLMs to use PDE descriptions and prior knowledge in design decisions\.

Lang\-PINN further adopts a multi\-agent framework that assigns PDE understanding, PINN architecture design, code generation, and execution feedback to different agents, while a Feedback Agent analyzes execution results and performs iterative corrections\([He et al\., 2026](https://arxiv.org/html/2609.23023#bib.bib10)\)\. These studies show that LLMs can not only understand PDE descriptions and propose PINN configurations, but can also use external execution results to further adjust generated outputs\.

However, real PINN training provides substantially richer information than the final performance score alone\. In addition to solution error, PDE residuals, boundary\- and initial\-condition errors, changes in different loss components, and optimization stability reveal the detailed behavior of a candidate configuration\. PINNsForge therefore focuses not merely on expanding the range of configurations that can be automated, but on how to continually incorporate such real training evidence into multi\-generation search so that the LLM can adapt subsequent designs according to successes and failures already observed on the current PDE\.

### 2\.3LLM\-Guided Evolutionary Algorithm Discovery

The role of LLMs in automated algorithm design has recently expanded from one\-shot program generation to iterative search\. FunSearch combines a pretrained LLM with an automated evaluator and evolutionary program search\. By repeatedly executing and selecting candidate programs, high\-quality programs are reused as context to guide subsequent generation\([Romera\-Paredes et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib30)\)\. This work demonstrates that real program execution results can provide an important signal for continual LLM\-driven design improvement\.

Evolution of Heuristics \(EoH\) further makes this idea explicit as an automated algorithm\-discovery framework that combines LLMs with evolutionary computation\([Liu et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib18);[Yao et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib51);[Liu et al\., 2026a](https://arxiv.org/html/2609.23023#bib.bib20)\)\. Its core idea is to embed the LLM in a “generate—execute—evaluate—evolve” loop, such that the model not only generates candidates but also continually proposes new algorithm designs based on high\-performing individuals and actual evaluation results\. ReEvo uses performance differences among candidates and LLM reflection to produce language\-based search feedback that guides subsequent heuristic improvement\([Ye et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib52)\), while LLaMEA uses runtime evaluation, selection, and iterative modification to continually generate and optimize candidate metaheuristics\([van Stein and Bäck, 2025](https://arxiv.org/html/2609.23023#bib.bib35);[van Stein et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib36)\)\.

Together, these methods suggest that, compared with random mutation or fixed numerical perturbations, an LLM can perform a more semantic form of evolutionary operation: it can understand the structure of existing candidates and their execution behavior, and use this information to decide which mechanisms should be inherited, modified, combined, or re\-explored\. This capability is particularly well suited to automated PINN design because executing a PINN candidate produces not only a scalar performance value but also training evidence with explicit physical meaning\.

PINNsForge further applies this “execute—feedback—evolve” principle to PINN configuration search\. Rather than determining the next search direction solely from final error, PINNsForge jointly uses PDE\-related priors, execution experience accumulated during the current search, and performance\-, physics\-, and optimization\-level diagnostics to provide the LLM with richer design evidence\. During evolution, different search strategies perform conservative modifications, physics\- and constraint\-guided corrections, network and optimization adjustments, multi\-parent synthesis, and novelty exploration, allowing configuration changes to coordinate around different problems exposed by actual training rather than being restricted to random perturbations of independent parameters\.

## 3Method

PINNsForge applies the “generate—execute—evaluate—evolve” principle of Evolution of Heuristics \(EoH\) to automated PINN design\. Rather than having an LLM produce a network architecture or training configuration only once, PINNsForge executes candidate configurations on the PDE and continually feeds performance, physical, and optimization evidence from real training back into subsequent search, so that new designs are grounded in observed solution behavior\.

Given a PDE, PINNsForge first combines equation characteristics with relevant PINN priors to generate an initial set of candidates with different design preferences\. After real PINN training, the system accumulates solution errors, PDE residuals, constraint satisfaction, and optimization behavior, and maintains the best\-performing Global Top\-3 over the entire search history\. The LLM then jointly modifies and recombines multiple coupled PINN design components based on high\-performing parents, execution feedback, and PDE knowledge through five complementary search strategies: Refine, Physics\-Guided, Architecture\-Guided, Synthesis, and Novelty\.

The overall search process can be summarized as

PDE Prior→Candidate Design→Real PINN Training→Execution Evidence→Joint Evolution\\boxed\{\\text\{PDE Prior\}\\rightarrow\\text\{Candidate Design\}\\rightarrow\\text\{Real PINN Training\}\\rightarrow\\text\{Execution Evidence\}\\rightarrow\\text\{Joint Evolution\}\}This closed loop continues until the search terminates\. The Global Top\-3 obtained from low\-fidelity search are then independently trained at high fidelity using new random initializations, and the best result is selected\. Figure 1 provides an overview of PINNsForge\.

![Refer to caption](https://arxiv.org/html/2609.23023v1/pinnsforge_framework.png)Figure 1:Overview of the PINNsForge framework\.### 3\.1Problem Definition

Letzzdenote the input coordinates of a PDE, and letu⁡\(z\)∈ℝduu\(z\)\\in\\mathbb\{R\}^\{d\_\{u\}\}denote the physical field to be solved over domain𝒟\\mathcal\{D\}, wheredu=1d\_\{u\}=1corresponds to a scalar field anddu\>1d\_\{u\}\>1to a vector field containing multiple state variables\. The target PDE can be written uniformly as

𝒩​\[u\]​\(z\)=0,z∈𝒟,\\mathcal\{N\}\[u\]\(z\)=0,\\qquad z\\in\\mathcal\{D\},
where𝒩\\mathcal\{N\}denotes the differential operator defined by the specific PDE\.

In addition to the governing equation, the problem typically includes initial conditions, boundary conditions, or other physical constraints, which can be written uniformly as

𝒞k​\[u\]​\(z\)=gk​\(z\),z∈Γk,\\mathcal\{C\}\_\{k\}\[u\]\(z\)=g\_\{k\}\(z\),\\qquad z\\in\\Gamma\_\{k\},
where𝒞k\\mathcal\{C\}\_\{k\}denotes thekk\-th constraint operator andΓk\\Gamma\_\{k\}denotes the corresponding region or boundary on which the constraint is imposed\.

A PINN uses a neural networkuθ​\(z\)u\_\{\\theta\}\(z\)to approximate the true physical fieldu⁡\(z\)u\(z\), and optimizes the network parametersθ\\thetaby minimizing a training objective composed of the PDE residual and physical constraints:

ℒ⁡\(θ\)=λr​ℒr\+∑kλk​ℒk,\\mathcal\{L\}\(\\theta\)=\\lambda\_\{r\}\\mathcal\{L\}\_\{r\}\+\\sum\_\{k\}\\lambda\_\{k\}\\mathcal\{L\}\_\{k\},
whereℒr\\mathcal\{L\}\_\{r\}denotes the PDE residual loss,ℒk\\mathcal\{L\}\_\{k\}denotes the loss term corresponding to an initial condition, boundary condition, or other physical constraint, andλr\\lambda\_\{r\}andλk\\lambda\_\{k\}are the associated weights\.

To evaluate the accuracy with which a PINN approximates the true physical field, we use mean squared error \(MSE\) over fixed evaluation pointszii=1Ne\{z\_\{i\}\}\_\{i=1\}^\{N\_\{e\}\}as the primary evaluation metric:

MSE=1Ne​∑i=1Ne‖uθ​\(zi\)−u⋆​\(zi\)‖22,\\operatorname\{MSE\}=\\frac\{1\}\{N\_\{e\}\}\\sum\_\{i=1\}^\{N\_\{e\}\}\\left\\\|u\_\{\\theta\}\(z\_\{i\}\)\-u^\{\\star\}\(z\_\{i\}\)\\right\\\|\_\{2\}^\{2\},
whereu⋆​\(zi\)u^\{\\star\}\(z\_\{i\}\)denotes the reference solution anduθ​\(zi\)u\_\{\\theta\}\(z\_\{i\}\)denotes the prediction of the trained PINN\. A lower MSE indicates a smaller discrepancy between the predicted and reference solutions\.

### 3\.2From PDE Priors to Real Execution Experience

PINNsForge jointly uses PDE\-related prior knowledge and candidate execution experience during search\. The former describes general physical principles and PINN design guidelines, providing directional guidance for candidate design; the latter comes from the actual training results of candidate configurations on the current PDE and is used to evaluate the effectiveness of different designs and analyze their performance bottlenecks\. As the search progresses, candidate generation gradually shifts from initial knowledge\-guided candidate generation toward execution\-grounded adaptation supported by real execution evidence\.

#### 3\.2\.1PDE\-Related Prior Knowledge

For a target PDE𝒫\\mathcal\{P\}, PINNsForge first analyzes its basic problem characteristics, including time dependence, periodicity, multiscale structure, geometric complexity, locally steep regions, and types of physical constraints, and retrieves relevant PINN design knowledge accordingly\.

This knowledge does not directly prescribe a specific PINN configuration\. Instead, it provides the relevant physical and numerical rationale to help the LLM understand the relationship between PDE characteristics and PINN design and formulate corresponding design hypotheses for the current problem\. Prior knowledge therefore does not restrict the search to a predefined set of configurations; rather, it is used throughout the search together with the accumulating execution evidence to guide candidate generation and adjustment in every round\.

#### 3\.2\.2Building Execution Experience from Real Training

Prior knowledge captures general design principles, while the real training results of candidates provide empirical evidence specific to the current PDE\. After each candidate finishes training, PINNsForge extracts execution evidence from three perspectives: performance evidence, including the candidate’s final solution accuracy and its performance relative to other candidates; physical evidence, including satisfaction of the PDE residual, initial conditions, boundary conditions, and other physical constraints; and optimization evidence, including loss reduction, training stagnation, gradient or loss imbalance, and numerical stability\. Together, these signals characterize the candidate’s actual training behavior and are used to analyze the main factors limiting its performance\.

PINNsForge organizes the execution evidence accumulated throughout the search into a run\-scoped posterior memory\. Unlike generally applicable PDE priors, this memory summarizes only the designs and execution outcomes actually observed for the current PDE within the current independent search run, providing problem\-specific empirical evidence for subsequent candidate generation and adjustment\.

Accordingly, as the search progresses, the information available to the LLM gradually expands to

\[PDE Prior\+Observed Execution Experience,\]

where PDE priors continue to provide general physical and design rationale, while execution experience continually adds evidence actually observed on the current PDE\. Subsequent candidate design can therefore follow general design principles while increasingly adapting to the actual training results obtained for the current PDE\.

### 3\.3LLM\-Driven Configuration Generation and Evolutionary Search

PINNsForge uses an LLM to generate complete PINN configurations rather than selecting candidates from a finite, pre\-enumerated configuration set\. For each generation, the LLM jointly receives relevant PDE prior knowledge, real training feedback from candidates, posterior experience accumulated during the run, and the current search role, and generates a new candidate configuration accordingly\. Candidates are represented as structured AlgorithmSpecs rather than arbitrary executable code; generated outputs are mapped to existing PINN implementations only after validity, compatibility, and budget checks\. We define an AlgorithmSpec as a structured representation of a complete PINN design that uniformly describes its representation, network architecture, sampling strategy, constraint handling, loss construction, optimization strategy, and training control\.

Design DimensionRepresentative Design ChoicesRepresentationCoordinate transforms, Fourier/periodic features, multiscale representationsArchitectureNetwork depth, width, activation functions, residual connectionsSamplingCollocation\-point allocation, sampling distributions, adaptive resamplingConstraintsSoft constraints, exact constraints, periodic or interface handlingLossPDE residual construction, loss weighting, normalization, causal weightingOptimizationOptimizer, learning rate, gradient handling, staged optimizationTraining ControlResampling, curriculum, checkpointing, and other training controls#### 3\.3\.1Knowledge\-Guided Initial Candidate Generation

At Generation 0, no candidate execution history exists for the current PDE, so the initial population is generated primarily from the PDE description and its relevant prior knowledge\.

To increase the design coverage of the initial population, PINNsForge defines eight initialization roles focusing on:

1. 1\.Robust reference design;
2. 2\.Physical constraints and temporal causality;
3. 3\.Local structure and adaptive sampling;
4. 4\.Network representation;
5. 5\.PDE residual and loss construction;
6. 6\.Optimization and numerical stability;
7. 7\.Model capacity and network topology;
8. 8\.Novel cross\-component combinations\.

These roles specify only the design focus during candidate generation and do not correspond to fixed algorithm templates\. Each role receives the same PDE description and prior knowledge and independently generates a complete configuration\. The resulting Generation 0 forms an initial population with diverse design preferences and serves as the starting point for subsequent evolutionary search\.

#### 3\.3\.2Candidate Execution and Feedback

Each generated candidate is evaluated through actual PINN training\. For candidate configurationaia\_\{i\}and target PDE𝒫\\mathcal\{P\}, we denote the execution process as

Execute⁡\(ai;𝒫\)→\(yi,Ei\),\\operatorname\{Execute\}\(a\_\{i\};\\mathcal\{P\}\)\\rightarrow\(y\_\{i\},E\_\{i\}\),
whereyiy\_\{i\}denotes the primary performance metric used for candidate comparison, andEiE\_\{i\}denotes the execution evidence obtained during training, including MSE, PDE residual, initial\- and boundary\-condition errors, loss evolution, optimization stagnation, and numerical stability\.

Candidate evaluation therefore produces not only a scalar performance metric for ranking, but also diagnostic information associated with training behavior\. The former is used for candidate selection, while the latter is used during subsequent candidate generation to identify which design components in the current configuration may need adjustment\.

#### 3\.3\.3Global Historical Elite Selection

Letℋg\\mathcal\{H\}\_\{g\}denote all candidates evaluated up to generationgg\. PINNsForge selects the three best distinct candidates from the entire search history rather than only from the previous generation:

𝒫g\+1=GlobalDistinctTop3⁡\(ℋg\)\.\\mathcal\{P\}\_\{g\+1\}=\\operatorname\{GlobalDistinctTop3\}\\left\(\\mathcal\{H\}\_\{g\}\\right\)\.
where𝒫g\+1\\mathcal\{P\}\_\{g\+1\}serves as the parent set for generating candidates in the next generation\. With global historical selection, high\-quality configurations found early in the search can continue to participate in subsequent evolution as long as they remain competitive, preventing them from being discarded solely through generational replacement\.

The Global Top\-3 simultaneously retains multiple validated high\-quality designs, providing a parent pool for subsequent local refinement and cross\-candidate synthesis\.

#### 3\.3\.4Execution\-Feedback\-Driven Joint Evolution

Starting from Generation 1, next\-generation candidates are generated from the Global Top\-3, their associated execution evidence, PDE priors, and different evolutionary strategies\.

Conventional parameter\-level mutation typically changes only a small number of independent variables, whereas representation, network architecture, sampling, constraints, loss design, and optimization are strongly coupled in PINNs\. For example, insufficient representation capacity may be jointly related to network capacity, loss balancing, and optimization, so adjusting a single hyperparameter may not improve training\. PINNsForge therefore allows the LLM to modify multiple related components simultaneously according to execution evidence\.

To cover different search directions, each generation uses five evolutionary strategies:

Search StrategyPrimary Evolution DirectionRefineUse the current best candidate as the primary parent, preserve effective existing designs, and make limited adjustments targeting its main shortcomingsPhysics\-GuidedAdjust sampling\-, constraint\-, and loss\-related designs based primarily on the PDE residual and initial/boundary\-condition errorsArchitecture\-GuidedAdjust network representation, architecture, and related training configurations based on representation capacity, model capacity, and optimization behaviorSynthesisCombine well\-performing design mechanisms from multiple candidates into a new unified configurationNoveltyExplore design combinations not yet covered by the current high\-performing candidates or existing search trajectoryThe strategies above define different search preferences\. Which parent designs are preserved, which components are modified, and the magnitude of each modification are determined by the LLM based on PDE characteristics, parent configurations, and execution evidence\.

The next\-generation candidate generation process can be written as

ag\+1=LLMGenerate⁡\(𝒫g\+1,E𝒫g\+1,K𝒫,s\),a\_\{g\+1\}=\\operatorname\{LLMGenerate\}\\left\(\\mathcal\{P\}\_\{g\+1\},E\_\{\\mathcal\{P\}\_\{g\+1\}\},K\_\{\\mathcal\{P\}\},s\\right\),
where𝒫g\+1\\mathcal\{P\}\_\{g\+1\}denotes the Global Top\-3 parents,E𝒫g\+1E\_\{\\mathcal\{P\}\_\{g\+1\}\}denotes the associated execution evidence,K𝒫K\_\{\\mathcal\{P\}\}denotes relevant prior knowledge for the target PDE, andssdenotes the evolutionary strategy used\.

Different strategies emphasize different forms of parent utilization and search directions\. Refine focuses on local improvement around the current best configuration; Physics\-Guided and Architecture\-Guided target issues related to physical constraints and model representation/optimization, respectively; Synthesis exploits complementary designs across multiple high\-performing candidates; and Novelty expands the coverage of the current search trajectory\. Through these mechanisms, PINNsForge extends candidate search from independent parameter adjustment to cross\-component configuration updates grounded in execution evidence\. New evidence produced by real training in each generation is added back to the posterior memory and, together with the global historical elites, forms the context for the next round of candidate generation, yielding a continual “generate—execute—evaluate—evolve” loop\.

### 3\.4Multi\-Fidelity Search

Execution\-feedback\-driven evolutionary search requires evaluating many candidates, while fully training every candidate would incur substantial computational cost\. Inspired by prior multi\-fidelity neural architecture search work, MF\-NAS\([Phan and Luong, 2024](https://arxiv.org/html/2609.23023#bib.bib25)\), we adopt a two\-stage strategy that combines low\-fidelity search with high\-fidelity evaluation\. Multi\-fidelity optimization typically allocates a small initial budget to many candidates and reserves larger computational budgets for a small number of promising candidates, balancing search efficiency and evaluation reliability\([Won et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib45)\)\.

At the same time, candidate rankings obtained under limited budgets do not necessarily reflect final performance after full training\. Prior NAS studies have shown that proxy evaluations used during search can produce rankings that are inconsistent with the true performance of independently trained candidates\([Yu et al\., 2020](https://arxiv.org/html/2609.23023#bib.bib57)\)\. PINNsForge therefore does not directly take the best low\-fidelity candidate as the final result, but instead divides the overall procedure into two stages:

Uniform 1K\-step Evolutionary Search⏟Low\-fidelity Search→Independent 10K\-step Evaluation of Global Top\-3⏟High\-fidelity Evaluation\\boxed\{\\underbrace\{\\text\{Uniform 1K\-step Evolutionary Search\}\}\_\{\\text\{Low\-fidelity Search\}\}\\;\\rightarrow\\;\\underbrace\{\\text\{Independent 10K\-step Evaluation of Global Top\-3\}\}\_\{\\text\{High\-fidelity Evaluation\}\}\}
The low\-fidelity stage rapidly advances evolutionary search under a uniform \(1\{,\}000\)\-step budget while continually updating the Global Top\-3\. In the high\-fidelity stage, each candidate in the final Global Top\-3 is reinitialized and independently trained for \(10\{,\}000\) optimizer steps\.

#### 3\.4\.1Low\-Fidelity Search

In the low\-fidelity search stage, every valid candidate is trained using a uniform budget of \(1\{,\}000\) optimizer steps\. For each candidate, the system records the PDE residual, boundary/initial\-condition errors, and convergence and stability information from optimization, and feeds these execution results into subsequent candidate generation\. The goal of the low\-fidelity stage is not to accurately predict a candidate’s final performance after full training, but to rapidly compare different PINN designs under a fixed computational budget, continually eliminate weaker designs, and retain more promising configurations for subsequent evolution\. Generation 0 contains 8 initial candidates, followed by at most 9 generations with 5 new candidates per generation\. Throughout the search, the system continually maintains the Global Top\-3 over the entire search history and uses them as an important reference for subsequent candidate generation\. As evolution proceeds, new candidates continually compete with the best historical configurations, enabling the search to explore and improve PINN designs under a limited budget\. When little improvement is observed for several consecutive generations, the system increases the intensity of novelty exploration; if further search still fails to produce effective improvement, low\-fidelity search terminates early and the final Global Top\-3 is passed to high\-fidelity evaluation\.

#### 3\.4\.2High\-Fidelity Evaluation

Low\-fidelity search continually explores and improves PINN designs over multiple generations while retaining the best\-performing Global Top\-3 from the entire search history\. The high\-fidelity stage fully trains these high\-quality candidates selected through evolution and determines the final best configuration for the current PDE\. For each selected candidate, the system keeps the PINN design identified during low\-fidelity search unchanged and reruns only the training process: 1\. Keep the candidate configuration unchanged and reconstruct the corresponding PINN; 2\. Use a new random initialization; 3\. Independently train from scratch for \(10\{,\}000\) optimizer steps; 4\. Select the candidate with the lowest MSE under high\-fidelity evaluation as the final output\.

## 4Experiments

In this section, we evaluate PINNsForge on a large and heterogeneous set of PDEs and focus on four questions:

1. 1\.Overall performance:Can PINNsForge automatically discover PINN configurations with higher solution accuracy across different types of PDEs?
2. 2\.Search cost:Can these performance gains be achieved under controlled search and training budgets?
3. 3\.Mechanistic contributions:What roles are played by PDE\-related knowledge, real execution feedback, and parent\-driven evolution?
4. 4\.Search process:Does the final performance genuinely arise from continual multi\-generation improvement driven by execution feedback, rather than from an accidentally strong configuration generated in Generation 0?

### 4\.1Experimental Setup

##### LLM backbone\.

Across all PINNsForge experiments, we use Qwen3\.6\-27B as the fixed LLM backbone for candidate PINN configuration generation, feedback analysis, and subsequent evolutionary search, without any additional PDE\-specific fine\-tuning\. Candidate generation uses a fixed role\-conditioned temperature schedule, with different sampling temperatures for different search roles; complete settings are provided in Appendix C\.2\.

##### Benchmark\.

The experiments include 25 PDE problems, 21 of which are drawn from the PINNacle benchmark\([Hao et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib9)\), spanning nonlinear transport, reaction–diffusion, multiscale and long\-time dynamics, complex geometries, high\-dimensional PDEs, inverse problems, and fluid systems\. Beyond the PINNacle problems, we additionally include four classical PDEs—Allen–Cahn, variable\-coefficient Darcy flow, the shallow\-water system, and Kovasznay flow\([Allen and Cahn, 1979](https://arxiv.org/html/2609.23023#bib.bib1);[Darcy, 1856](https://arxiv.org/html/2609.23023#bib.bib3);[de Saint\-Venant, 1871](https://arxiv.org/html/2609.23023#bib.bib32);[Kovasznay, 1948](https://arxiv.org/html/2609.23023#bib.bib15)\)—to broaden the range of equation structures\. Definitions of all 25 PDEs are provided in Appendix A\.

##### Baselines\.

We select RoPINN, PINNsFormer, and PINNsAgent\([Wu et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib47);[Zhao et al\., 2024](https://arxiv.org/html/2609.23023#bib.bib61);[Wuwu et al\., 2025](https://arxiv.org/html/2609.23023#bib.bib49)\)as the unified comparison baselines, representing PINN optimization, architectural design, and LLM\-driven automated configuration, respectively\. We additionally report PINNacle Best as an oracle\-style reference\. Because it combines the best method across different tasks, it is excluded from average\-rank calculations and statistical significance tests\.

##### Evaluation Metrics\.

For each PDE, before the experiments begin, we fix a partition of the evaluation points associated with the reference solution into two non\-overlapping sets: 30% as the search\-validation reference set and 70% as the final\-test reference set\. The search\-validation reference set is used only for candidate comparison, Global Top\-3 updates, and stagnation and escape decisions during PINNsForge search\. The final\-test reference set remains unseen throughout search and model selection and is used only for performance evaluation after the final candidate has been determined\. We use solution MSE on the final\-test set as the primary reported metric and report the mean±\\pmstandard deviation over 10 independent runs\.

##### Training Protocol\.

RoPINN, PINNsFormer, and PINNsAgent are trained according to their respective unified experimental settings, and final performance is evaluated on the same final\-test reference set\. PINNsForge uses multi\-fidelity search\. During the low\-fidelity stage, each candidate is trained for1​K1\\mathrm\{K\}optimizer steps, and the Global Top\-3 over the full search history are continually maintained according to search\-validation MSE\. After search, the Top\-3 candidate configurations are kept unchanged and each is trained for10​K10\\mathrm\{K\}optimizer steps from a new random initialization to obtain the final result\. The complete PINNsForge search protocol is provided in Appendix C\.1\.

##### Implementation Details\.

All experiments are run on a Linux server with an NVIDIA GeForce RTX 4090 GPU, Python 3\.12, PyTorch 2\.6\.0\+cu124, and CUDA 12\.4\. Each PDE search uses a single GPU and FP32 training\. Random seeds control network initialization and training sampling; complete evaluation and randomness settings are provided in Appendix C\.3\.

### 4\.2Main Results

Table 1 summarizes the main experimental results on 25 PDE problems\. All results are obtained from 10 independent runs and report the mean±\\pmstandard deviation of the final selected result\. Boldface denotes the lowest mean MSE among the four unified comparison methods: RoPINN, PINNsFormer, PINNsAgent, and PINNsForge\.Rankdenotes the per\-task rank of PINNsForge among these four methods\. PINNacle Best reports the best result provided by the PINNacle benchmark for each PDE and is included only as a reference; it is excluded from rankings and statistical tests among the unified methods\.

Table 1: MSE on 25 PDE problems\.

PDE ProblemRoPINNPINNsFormerPINNsAgentPINNsForgeRankPINNacle Bestburgers\_1d7\.31E\-04±\\pm6\.13E\-051\.08E\-03±\\pm8\.21E\-066\.52E\-05±\\pm1\.55E\-056\.50E\-05±\\pm4\.42E\-0717\.90E\-05burgers\_2d2\.40E\-01±\\pm3\.12E\-022\.65E\-01±\\pm6\.83E\-021\.98E\-01±\\pm1\.66E\-021\.95E\-01±\\pm1\.31E\-0211\.09E\-01gray\_scott\_2d8\.64E\-03±\\pm4\.53E\-051\.57E\-02±\\pm1\.32E\-034\.36E\-03±\\pm3\.08E\-053\.60E\-03±\\pm2\.65E\-0414\.32E\-03heat\_2d\_complex\_geometry4\.00E\-04±\\pm9\.10E\-052\.78E\-03±\\pm5\.77E\-051\.84E\-03±\\pm9\.95E\-044\.45E\-04±\\pm9\.72E\-0528\.53E\-04heat\_2d\_long\_time2\.56E\+00±\\pm4\.53E\-013\.83E\+00±\\pm1\.64E\-014\.32E\+00±\\pm2\.78E\+001\.13E\+00±\\pm4\.21E\-0211\.13E\+00heat\_2d\_multiscale2\.75E\-05±\\pm7\.52E\-065\.57E\-05±\\pm9\.93E\-063\.51E\-05±\\pm2\.20E\-051\.19E\-05±\\pm2\.16E\-0615\.27E\-05heat\_2d\_varying\_coefficient6\.24E\-01±\\pm4\.53E\-031\.19E\-02±\\pm3\.42E\-035\.34E\-03±\\pm3\.85E\-031\.59E\-03±\\pm4\.77E\-0411\.76E\-03heat\_5d5\.39E\-07±\\pm1\.31E\-071\.87E\-07±\\pm1\.91E\-083\.34E\-07±\\pm7\.65E\-083\.21E\-08±\\pm3\.26E\-0818\.52E\+00heat\_inverse\_2d4\.31E\-02±\\pm7\.31E\-038\.89E\-02±\\pm4\.22E\-026\.32E\-02±\\pm8\.63E\-027\.15E\-03±\\pm4\.91E\-0315\.66E\-03kuramoto\_sivashinsky\_1d6\.54E\+00±\\pm5\.14E\-015\.13E\+00±\\pm4\.51E\-031\.09E\+00±\\pm3\.18E\-021\.03E\+00±\\pm8\.40E\-0211\.04E\+00navier\_stokes\_2d\_backstep9\.64E\-03±\\pm5\.13E\-036\.02E\-04±\\pm4\.33E\-054\.96E\-03±\\pm1\.45E\-032\.63E\-04±\\pm4\.71E\-0513\.37E\-04navier\_stokes\_2d\_classic2\.88E\-05±\\pm3\.61E\-061\.49E\-05±\\pm1\.26E\-058\.48E\-06±\\pm6\.90E\-062\.48E\-06±\\pm2\.31E\-0612\.33E\-05navier\_stokes\_2d\_long\_time5\.43E\+02±\\pm5\.63E\+007\.43E\+02±\\pm1\.31E\+025\.75E\+02±\\pm7\.63E\+005\.04E\+02±\\pm4\.95E\+0015\.05E\+02poisson\_2d\_classic8\.57E\-02±\\pm6\.74E\-022\.17E\-01±\\pm7\.65E\-026\.34E\-01±\\pm1\.43E\-011\.61E\-02±\\pm3\.95E\-0415\.00E\-05poisson\_boltzmann\_2d7\.91E\-01±\\pm7\.14E\-011\.34E\-02±\\pm4\.31E\-036\.97E\-02±\\pm7\.43E\-025\.38E\-03±\\pm9\.13E\-0416\.99E\-05poisson\_2d\_many\_area2\.89E\+00±\\pm3\.51E\-015\.36E\+00±\\pm3\.79E\+003\.00E\+00±\\pm1\.01E\+002\.21E\+00±\\pm3\.56E\-0111\.83E\+00poisson\_3d\_complex\_geometry2\.16E\-01±\\pm1\.53E\-028\.16E\-03±\\pm5\.58E\-041\.60E\-02±\\pm1\.08E\-021\.04E\-03±\\pm2\.34E\-0419\.51E\-04poisson\_5d4\.09E\-07±\\pm7\.43E\-082\.92E\-07±\\pm1\.99E\-072\.01E\-06±\\pm1\.09E\-065\.64E\-08±\\pm3\.95E\-0812\.09E\-06poisson\_inverse\_2d5\.32E\-03±\\pm3\.46E\-031\.51E\-03±\\pm8\.65E\-042\.41E\-02±\\pm4\.73E\-022\.57E\-04±\\pm1\.29E\-0411\.23E\-04wave\_1d4\.19E\-02±\\pm4\.38E\-033\.61E\-02±\\pm6\.65E\-033\.25E\-02±\\pm3\.45E\-023\.18E\-02±\\pm5\.73E\-0313\.01E\-03wave\_2d\_heterogeneous4\.63E\-02±\\pm6\.73E\-038\.37E\-02±\\pm1\.18E\-025\.17E\-02±\\pm8\.03E\-032\.84E\-02±\\pm7\.69E\-0412\.99E\-02allen\_cahn\_1d7\.45E\-08±\\pm4\.61E\-086\.02E\-09±\\pm1\.08E\-091\.44E\-08±\\pm3\.50E\-091\.13E\-09±\\pm2\.26E\-101–darcy\_flow\_2d3\.48E\-09±\\pm1\.63E\-091\.35E\-10±\\pm2\.90E\-115\.33E\-07±\\pm2\.65E\-071\.64E\-11±\\pm3\.28E\-121–shallow\_water\_2d5\.78E\-07±\\pm7\.13E\-072\.24E\-07±\\pm8\.44E\-086\.81E\-07±\\pm1\.25E\-077\.15E\-08±\\pm2\.86E\-081–kovasznay\_flow\_2d6\.13E\-07±\\pm4\.31E\-073\.80E\-08±\\pm9\.63E\-093\.82E\-07±\\pm3\.66E\-071\.77E\-09±\\pm4\.43E\-101–Friedman Avg\. Rank3\.042\.962\.961\.04––Wilcoxon vs\. PINNsForge\+\+\+–––\+,=, and\-indicate that PINNsForge is significantly better than, statistically indistinguishable from, or significantly worse than the corresponding method according to the Wilcoxon signed\-rank test with Holm correction atα=0\.05\\alpha=0\.05\.

#### 4\.2\.1Overall Performance

PINNsForge achieves the lowest mean MSE among the four unified comparison methods on24 of 25PDEs, ranking second only onheat\_2d\_complex\_geometry\. Using the best\-performing baseline on each PDE as the reference, PINNsForge yields amedian MSE reduction of 59\.85%, while the geometric mean of the best\-baseline\-to\-PINNsForge MSE ratio reaches2\.78×\\times\. Its overall advantage therefore does not arise from extreme gains on a small number of tasks, but from relatively consistent improvements across tasks\.

These gains span PDEs with diverse structures and solution characteristics\. PINNsForge reduces MSE by95\.34%, 87\.85%, and 87\.25%on Kovasznay flow, variable\-coefficient Darcy flow, and 3D complex\-geometry Poisson, respectively, and by83\.41%, 82\.98%, and 82\.83%on inverse heat, inverse Poisson, and 5D heat\. Moreover, PINNsForge achieves the lowest mean MSE among the unified methods on all four additional PDEs outside PINNacle\. These results indicate that the performance gains are not limited to a particular equation family, but appear across fluid systems, variable\-coefficient problems, complex geometries, inverse problems, and high\-dimensional PDEs\.

#### 4\.2\.2Stability and Statistical Significance

In terms of relative variation across repeated runs, PINNsForge achieves lower MSE without an obvious loss of run\-to\-run stability\. Measured by the coefficient of variation \(CV; standard deviation / mean\), PINNsForge has a median CV of0\.182across the 25 problems, comparable to0\.177for RoPINN and0\.184for PINNsFormer, and notably lower than0\.497for PINNsAgent\.

We further apply a Friedman test to the overall rankings of the four unified methods across the 25 PDEs, which indicates a significant difference among the methods:

χF2=42\.696,p=2\.86×10−9\.\\chi\_\{F\}^\{2\}=42\.696,\\qquad p=2\.86\\times 10^\{\-9\}\.
PINNsForge has an average rank of1\.04, compared with3\.04, 2\.96, and 2\.96for RoPINN, PINNsFormer, and PINNsAgent, respectively\. Given the significant Friedman test, we further conduct paired Wilcoxon signed\-rank tests onlog10⁡\(MSE\)\\log\_\{10\}\(\\mathrm\{MSE\}\)across the 25 PDEs and apply Holm correction for multiple comparisons\. PINNsForge shows a statistically significant advantage over RoPINN, PINNsFormer, and PINNsAgent, with the adjusted significance level for all three comparisons given by

padj=1\.79×10−7\.p\_\{\\mathrm\{adj\}\}=1\.79\\times 10^\{\-7\}\.
Taken together, the per\-task win rate, relative error reduction, stability across repeated runs, and cross\-task statistical tests show that the performance gains of PINNsForge are broad and consistent across PDEs rather than being driven by a small number of problems\.

### 4\.3Component Ablation

To analyze the contribution of different mechanisms in PINNsForge to final performance, we construct three ablated variants on 12 representative PDEs\.

- •w/o Knowledge Base: Remove the PDE/PINN prior knowledge provided to the LLM and the posterior knowledge formed during the run, to evaluate the role of knowledge guidance in candidate generation and subsequent search\.
- •w/o Execution Feedback: Retain PDE prior knowledge and parent configurations, but do not provide subsequent LLM calls with the current candidates’ PDE residuals, constraint errors, optimization diagnostics, or run\-scoped posterior, to evaluate the contribution of real training\-result feedback to subsequent design\.
- •w/o Evolutionary Search: Remove parent\-based evolution so that subsequent candidates are generated independently and do not inherit, modify, or recombine existing configurations from the Global Top\-3\. This ablation tests whether the performance gains arise from continual evolution rather than simply evaluating more candidates under the same budget\.

Table 2: Component ablation settings for PINNsForge\.

VariantPrior KnowledgeRun\-scoped PosteriorExecution Feedback to LLMParent\-based Evolutionw/o Knowledge Base×\\times×\\times✓\\checkmark✓\\checkmarkw/o Execution Feedback✓\\checkmark×\\times×\\times✓\\checkmarkw/o Evolutionary Search✓\\checkmark✓\\checkmark✓\\checkmark×\\timesPINNsForge✓\\checkmark✓\\checkmark✓\\checkmark✓\\checkmarkBecause MSE values span multiple orders of magnitude across PDEs, we additionally compute the error ratio of each ablated variant relative to the full PINNsForge on each task:

r=MSEablationMSEfull,r=\\frac\{\\mathrm\{MSE\}\_\{\\mathrm\{ablation\}\}\}\{\\mathrm\{MSE\}\_\{\\mathrm\{full\}\}\},
and use the medianrmedianr\_\{\\mathrm\{median\}\}across the 12 tasks to summarize the overall performance change for each ablated variant\.

Table 3: Component ablation results\. Lower MSE is better; the last row reports the median cross\-task error ratio relative to the full PINNsForge\.

PDE ProblemRepresentative Characteristicw/o Knowledge Basew/o Execution Feedbackw/o Evolutionary SearchPINNsForgeburgers\_1dNonlinear / shock7\.42E\-05±\\pm6\.31E\-056\.32E\-04±\\pm7\.82E\-052\.78E\-04±\\pm3\.24E\-056\.50E\-05±\\pm4\.42E\-07burgers\_2d2D nonlinear transport2\.44E\-01±\\pm2\.86E\-026\.71E\-01±\\pm8\.37E\-021\.67E\+00±\\pm2\.31E\-011\.95E\-01±\\pm1\.31E\-02gray\_scott\_2dCoupled reaction–diffusion5\.29E\-03±\\pm4\.91E\-045\.36E\-03±\\pm5\.62E\-036\.23E\-03±\\pm6\.87E\-043\.60E\-03±\\pm2\.65E\-04heat\_2d\_multiscaleMultiscale dynamics3\.81E\-05±\\pm4\.28E\-062\.08E\-04±\\pm2\.74E\-041\.30E\-03±\\pm1\.66E\-031\.19E\-05±\\pm2\.16E\-06heat\_2d\_varying\_coefficientVariable coefficient2\.18E\-03±\\pm2\.03E\-032\.57E\-03±\\pm2\.71E\-044\.15E\-03±\\pm4\.86E\-041\.59E\-03±\\pm4\.77E\-04heat\_2d\_complex\_geometryComplex geometry1\.90E\-03±\\pm2\.41E\-031\.59E\+00±\\pm2\.14E\-011\.97E\-02±\\pm2\.77E\-034\.45E\-04±\\pm9\.72E\-05poisson\_3d\_complex\_geometry3D complex geometry8\.79E\-03±\\pm9\.38E\-042\.00E\-02±\\pm2\.66E\-032\.49E\-02±\\pm3\.18E\-031\.04E\-03±\\pm2\.34E\-04poisson\_5dHigh\-dimensional PDE2\.55E\-07±\\pm2\.87E\-071\.24E\-07±\\pm1\.41E\-082\.06E\-05±\\pm2\.72E\-065\.64E\-08±\\pm3\.95E\-08poisson\_inverse\_2dInverse problem6\.94E\-04±\\pm7\.64E\-053\.72E\-03±\\pm4\.58E\-043\.14E\-04±\\pm3\.63E\-042\.57E\-04±\\pm1\.29E\-04navier\_stokes\_2d\_backstepNonlinear fluid dynamics2\.05E\-03±\\pm2\.68E\-044\.50E\-03±\\pm6\.13E\-033\.06E\-03±\\pm3\.97E\-032\.63E\-04±\\pm4\.71E\-05darcy\_flow\_2dHeterogeneous2\.76E\-10±\\pm3\.17E\-102\.97E\-09±\\pm4\.06E\-101\.07E\-08±\\pm1\.48E\-091\.64E\-11±\\pm3\.28E\-12shallow\_water\_2dHyperbolic system3\.18E\-07±\\pm3\.52E\-082\.32E\-07±\\pm2\.69E\-085\.59E\-07±\\pm7\.14E\-087\.15E\-08±\\pm2\.86E\-08Median error ratiormedianr\_\{\\mathrm\{median\}\}–3\.74×\\times12\.10×\\times10\.10×\\times1\.00×\\timesOverall, removing the Knowledge Base, Execution Feedback, and Evolutionary Search increases the median MSE across the 12 tasks to3\.74×\\times,12\.10×\\times, and10\.10×\\timesthat of the full PINNsForge, respectively\. The largest degradations occur when removing Execution Feedback and Evolutionary Search, indicating that feedback from real training and continual evolution from high\-performing parents are two key mechanisms underlying the stable performance gains of PINNsForge\.

The three ablations also exhibit distinct functional patterns\. The Knowledge Base has a relatively small effect on problems such as Burgers, but a larger effect on complex\-geometry, high\-dimensional PDE, backstep\-flow, and Darcy\-flow tasks\. This suggests that PDE priors primarily provide a better\-grounded starting point for search and reduce ineffective exploration\. Without Execution Feedback, errors increase substantially on most tasks; forheat\_2d\_complex\_geometry, the MSE rises from4\.45×10−44\.45\\times 10^\{\-4\}to \(1\.59\)\. This indicates that static knowledge available before search is insufficient for continual candidate improvement, whereas performance, physical, and optimization information produced by real training provides direct evidence specific to the current PDE for subsequent design\.

When Evolutionary Search is removed, the same candidate\-evaluation budget and the available knowledge and execution feedback are retained, but existing designs are no longer inherited, modified, or recombined from the Global Top\-3\. The median cross\-task MSE still increases to10\.10×\\timesthat of the full PINNsForge\. The degradation is particularly pronounced on multiscale, high\-dimensional, and complex\-geometry problems, indicating that preserving strong historical candidates and continuing to modify and recombine them is an important mechanism for translating accumulated execution experience into subsequent performance gains\.

Overall, the Knowledge Base provides a grounded starting point for search, Execution Feedback feeds real training behavior into subsequent design, and Evolutionary Search enables effective designs to be continually retained, modified, and recombined\. Together, these components form the execution\-feedback\-driven evolutionary search loop of PINNsForge\.

### 4\.4Search Process Analysis

To further analyze the search process itself, we select four representative PDEs and examine their multi\-generation search trajectories together with the performance changes between the low\- and high\-fidelity stages\. These problems cover nonlinear transport, multiscale dynamics, high\-dimensional PDEs, and fluid systems, allowing us to observe the search behavior of PINNsForge under different solution characteristics\. Search trajectories and corresponding results for the remaining PDEs are provided in Appendix D\.

![Refer to caption](https://arxiv.org/html/2609.23023v1/search_trajectories.png)Figure 2:Evolutionary Search Trajectories across Representative PDE Problems\.Figure 2 shows the evolutionary search trajectories for the four representative PDEs\. Each curve records the best low\-fidelity MSE found up to the current generation during search, while the star marker denotes the best result obtained after independently training the Global Top\-3 at high fidelity from new random initializations\.

Figure 2: Evolutionary Search Trajectories across Representative PDE Problems

Table 4 further summarizes one independent search trajectory for each of the four PDEs\. Generation\-0 Best denotes the lowest MSE among the initial candidates, Best LF denotes the best result found over the entire low\-fidelity search stage, and Final HF denotes the best MSE obtained after independently training the Global Top\-3 under the high\-fidelity budget\.

PDEGeneration\-0 Best MSEBest LF MSEFinal HF MSEG0→\\rightarrowBest LFG0→\\rightarrowFinal HFburgers\_1d3\.64E\-032\.05E\-04 \(G8\)6\.51E\-0517\.74×\\times55\.82×\\timesheat\_2d\_multiscale3\.01E\-029\.140E\-05 \(G9\)1\.014E\-05330\.08×\\times2,975\.09×\\timespoisson\_5d6\.44E\-061\.01E\-07 \(G9\)5\.65E\-0863\.58×\\times114\.03×\\timesnavier\_stokes\_2d\_classic1\.86E\-031\.61E\-05 \(G9\)2\.35E\-06115\.60×\\times792\.09×\\timesThe search trajectories show that all four representative problems discover a new global\-best candidate after Generation 0\. In particular,heat\_2d\_multiscale,poisson\_5d, andnavier\_stokes\_2d\_classiccontinue to produce new best configurations through Generation 9\. This indicates that final performance is not primarily determined by a one\-shot LLM generation during initialization; subsequent multi\-generation search continues to discover more promising PINN designs\.

Relative to the best Generation\-0 candidate, the best configurations eventually found during low\-fidelity search improve MSE by17\.74×17\.74\\times,330\.08×330\.08\\times,63\.58×63\.58\\times, and115\.60×115\.60\\times, respectively\. This shows that continual search combining strong historical candidates with real execution feedback can progressively discover more effective PINN configurations across different types of PDEs\.

Building on these results, the Global Top\-3 candidate configurations are kept unchanged and each is trained at high fidelity from a new random initialization\. The high\-fidelity stage does not generate new candidates; instead, it evaluates how the strong configurations selected during low\-fidelity search perform under more complete training and re\-establishes the final ranking\. Relative to Generation 0, the four representative runs ultimately achieve MSE improvements of55\.82×55\.82\\times,2,975\.09×2\{,\}975\.09\\times,114\.03×114\.03\\times, and792\.09×792\.09\\times, respectively\. These results indicate that low\-fidelity search effectively identifies promising configurations, while high\-fidelity training further confirms their performance after more complete optimization\.

### 4\.5Training Budget

To compare the training cost required by each method to obtain its final result, we count the optimizer\-step budget used for candidate selection on a single PDE\.

For RoPINN, PINNsFormer, and PINNsAgent, each independent experiment trains five candidates with different random seeds, using 20K optimizer steps per candidate, and selects the lowest MSE as the final result\. The total training budget per independent experiment is therefore

5×20​K=100​K5\\times 20K=100K
optimizer steps\.

PINNsForge uses multi\-fidelity budget allocation\. Generation 0 contains 8 candidates, followed by at most 9 generations with 5 candidates generated per generation, resulting in at most

low\-fidelity candidates\. Each candidate is trained for 1K optimizer steps, giving a maximum low\-fidelity budget of 53K\. After search, the Top\-3 candidates are selected from the full search history, their configurations are kept unchanged, and each is independently trained from a new random initialization for 10K optimizer steps for high\-fidelity evaluation, adding a further 30K steps\. The maximum total training budget of PINNsForge is therefore

53​K\+30​K=83​K53K\+30K=83K
optimizer steps, which is lower than the 100K budget used by the three unified comparison baselines\.

Table 5: Optimizer\-step budget for obtaining the final result on a single PDE\.

MethodResult Selection ProtocolOptimizer\-Step BudgetRoPINNBest of 5 independent 20K runs100KPINNsFormerBest of 5 independent 20K runs100KPINNsAgentBest of 5 independent 20K runs100KPINNsForge53×\\times1K LF \+ Top\-3×\\times10K HF≤83​K\\leq 83\\mathrm\{K\}

## 5Conclusion

We propose PINNsForge, an execution\-feedback\-driven automated search framework over the complete PINN design space\. The method jointly incorporates PDE prior knowledge, performance and physical diagnostics produced by real training, run\-scoped posterior experience, and high\-performing candidates from the full search history into a multi\-generation search process\. It uses an LLM to jointly modify and recombine coupled design components spanning network representation, architecture, sampling, constraints, loss design, and optimization, forming a closed\-loop “design—execute—feedback—evolve” search mechanism\. The experimental results show that automated PINN design can be extended from one\-shot configuration generation to cross\-component evolutionary search continually driven by real solution evidence, providing a unified and scalable implementation paradigm for automated PINN design across heterogeneous PDEs\. Across 25 heterogeneous PDEs, PINNsForge achieves the lowest mean MSE among the four unified comparison methods on 24 tasks, with an average rank of 1\.04; ablations on 12 representative problems further indicate that execution feedback and parent\-driven evolutionary search play key roles in overall performance\.

## 6Broader Impact

PINNsForge is intended to reduce the reliance of PINN design in scientific machine learning on manual expertise and repeated trial and error\. By preserving structured candidate configurations together with their corresponding execution evidence, it also improves the traceability and auditability of the automated solver\-design process\. The framework may support automated model design and method exploration across different types of PDEs\.

It is important to note that a low reference error obtained through automated search does not by itself guarantee physical correctness or generalization\. Candidate solvers may still degrade in unsampled regions, out\-of\-distribution regimes, or under more complex parameter settings, and may inadequately satisfy local physical constraints or conservation relations\. PINNsForge is therefore better viewed as a tool for assisting scientific computing and algorithm design, and its outputs should still be assessed through independent validation, physical\-consistency checks, uncertainty analysis, and domain\-expert review\.

In addition, LLM\-driven multi\-round search introduces extra inference and computational overhead\. Evaluations of such automated design methods should therefore consider not only final solution accuracy, but also training cost, LLM\-call overhead, and overall resource consumption to support fairer and more comprehensive comparisons\.

## References

- Allen and Cahn \(1979\)S\. M\. Allen and J\. W\. CahnA microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening\.Note:*Acta Metallurgica*, 27\(6\):1085–1095Cited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px2.p1.1)\.
- Baiet al\.\(2025\)J\. Bai, Y\. Wang, H\. Jeong, S\. Chu, Q\. Wang, L\. Alzubaidi, X\. Zhuang, T\. Rabczuk, Y\. M\. Xie, X\. Feng, and Y\. GuTowards the future of physics\- and data\-guided AI frameworks in computational mechanics\.Note:*Acta Mechanica Sinica*, 41\(7\):225340doi:10\.1007/s10409\-025\-25340\-xCited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- Caoet al\.\(2025\)F\. Cao, X\. Guo, X\. Dong, and D\. YuanwbPINN: Weight balanced physics\-informed neural networks for multi\-objective learning\.Note:*Applied Soft Computing*, 170:112632doi:10\.1016/j\.asoc\.2024\.112632Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Carrilloet al\.\(2024\)H\. Carrillo, T\. de Wolff, L\. Martí, and N\. Sanchez\-PiEvolutionary multi\-objective physics\-informed neural networks: The MOPINNs approach\.Note:*AI Communications*, 37\(3\):397–409doi:10\.3233/AIC\-230073Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p2.1)\.
- Darcy \(1856\)H\. DarcyLes fontaines publiques de la ville de Dijon: Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau\.Note:Victor Dalmont, ParisCited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px2.p1.1)\.
- de Saint\-Venant \(1871\)A\. J\. C\. B\. de Saint\-VenantThéorie du mouvement non\-permanent des eaux, avec application aux crues des rivières et à l’introduction des marées dans leur lit\.Note:*Comptes Rendus des Séances de l’Académie des Sciences*, 73:147–154, 237–240Cited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px2.p1.1)\.
- Franciset al\.\(2026\)J\. M\. Francis, C\. Godavarma, and N\. ChamakuriAn adaptive residual\-structured and scheduled loss weighting framework for enhancing physics\-informed neural networks in solving partial differential equations\.Note:*Communications in Nonlinear Science and Numerical Simulation*, 153:109531doi:10\.1016/j\.cnsns\.2025\.109531Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- Franciset al\.\(2025\)J\. M\. Francis, R\. Kumar, C\. G, and N\. ChamakuriS\-PINN: Physics\-informed neural networks for solving weakly hyperbolic systems\.Note:*Physica Scripta*, 100\(12\):126002doi:10\.1088/1402\-4896/ae2391Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- Gaoet al\.\(2025\)B\. Gao, R\. Yao, and Y\. LiPhysics\-informed neural networks with adaptive loss weighting algorithm for solving partial differential equations\.Note:*Computers & Mathematics with Applications*, 181:216–227doi:10\.1016/j\.camwa\.2025\.01\.007Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Haoet al\.\(2024\)Z\. Hao, J\. Yao, C\. Su, H\. Su, Z\. Wang, F\. Lu, Z\. Xia, Y\. Zhang, S\. Liu, L\. Lu, and J\. ZhuPINNacle: A comprehensive benchmark of physics\-informed neural networks for solving PDEs\.Note:*Advances in Neural Information Processing Systems*, 37Cited by:[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p1.1),[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px2.p1.1)\.
- Heet al\.\(2026\)X\. He, L\. You, H\. Tian, B\. Han, I\. Tsang, and Y\. OngLang\-PINN: From language to physics\-informed neural networks via a multi\-agent framework\.Note:*ICLR 2026 Workshop on AI with Recursive Self\-Improvement*Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p4.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p5.1)\.
- Jagtap and Karniadakis \(2020\)A\. D\. Jagtap and G\. E\. KarniadakisExtended physics\-informed neural networks \(XPINNs\): A generalized space\-time domain decomposition based deep learning framework for nonlinear partial differential equations\.Note:*Communications in Computational Physics*, 28\(5\):2002–2041Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Jagtapet al\.\(2020\)A\. D\. Jagtap, K\. Kawaguchi, and G\. E\. KarniadakisAdaptive activation functions accelerate convergence in deep and physics\-informed neural networks\.Note:*Journal of Computational Physics*, 404:109136Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Kaplarević\-Mališićet al\.\(2023\)A\. Kaplarević\-Mališić, B\. Andrijević, F\. Bojović, S\. Nikolić, L\. Krstić, B\. Stojanović, and M\. IvanovićIdentifying optimal architectures of physics\-informed neural networks by evolutionary strategy\.Note:*Applied Soft Computing*, 146:110646Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p2.1)\.
- Kharazmiet al\.\(2021\)E\. Kharazmi, Z\. Zhang, and G\. E\. Karniadakishp\-VPINNs: Variational physics\-informed neural networks with domain decomposition\.Note:*Computer Methods in Applied Mechanics and Engineering*, 374:113547Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Kovasznay \(1948\)L\. I\. G\. KovasznayLaminar flow behind a two\-dimensional grid\.Note:*Proceedings of the Cambridge Philosophical Society*, 44\(1\):58–62Cited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px2.p1.1)\.
- Krishnapriyanet al\.\(2021\)A\. S\. Krishnapriyan, A\. Gholami, S\. Zhe, R\. M\. Kirby, and M\. W\. MahoneyCharacterizing possible failure modes in physics\-informed neural networks\.Note:*Advances in Neural Information Processing Systems*, 34:26548–26560Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p4.1)\.
- Lin and Chen \(2025\)S\. Lin and Y\. ChenCausality\-guided adaptive sampling method for physics\-informed neural networks solving forward problems of partial differential equations\.Note:*Physica D: Nonlinear Phenomena*, 481:134878doi:10\.1016/j\.physd\.2025\.134878Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Liuet al\.\(2026a\)F\. Liu, Y\. Liu, Q\. Zhang, X\. Tong, and M\. YuanEoH\-S: Evolution of heuristic set using LLMs for automated heuristic design\.Note:*Proceedings of the AAAI Conference on Artificial Intelligence*, 40\(43\):37090–37098doi:10\.1609/aaai\.v40i43\.41038Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- Liuet al\.\(2024\)F\. Liu, X\. Tong, M\. Yuan, X\. Lin, F\. Luo, Z\. Wang, Z\. Lu, and Q\. ZhangEvolution of Heuristics: Towards efficient automatic algorithm design using large language model\.Note:In*Proceedings of the 41st International Conference on Machine Learning*, PMLR 235:32201–32223Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- Liuet al\.\(2026b\)F\. Liu, Y\. Yao, P\. Guo, Z\. Yang, X\. Lin, Z\. Zhao, X\. Tong, K\. Mao, Z\. Lu, Z\. Wang, M\. Yuan, and Q\. ZhangA systematic survey on large language models for algorithm design\.Note:*ACM Computing Surveys*, 58\(8\), Article 218:1–32doi:10\.1145/3787585Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Liuet al\.\(2025\)Z\. Liu, Y\. Liu, X\. Yan, W\. Liu, H\. Nie, S\. Guo, and C\. ZhangAutomatic network structure discovery of physics informed neural networks via knowledge distillation\.Note:*Nature Communications*, 16:9558doi:10\.1038/s41467\-025\-64624\-3Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1),[§1](https://arxiv.org/html/2609.23023#S1.p4.1)\.
- Luoet al\.\(2025\)J\. Luo, Y\. Yang, Y\. Yuan, S\. Xu, and W\. HaoAn imbalanced learning\-based sampling method for physics\-informed neural networks\.Note:*Journal of Computational Physics*, 534:114010doi:10\.1016/j\.jcp\.2025\.114010Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- McClenny and Braga\-Neto \(2023\)L\. D\. McClenny and U\. M\. Braga\-NetoSelf\-adaptive physics\-informed neural networks\.Note:*Journal of Computational Physics*, 474:111722Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Michaloglouet al\.\(2026\)A\. Michaloglou, I\. Papadimitriou, I\. Gialampoukidis, S\. Vrochidis, and I\. KompatsiarisPhysics\-Informed Neural Networks in Materials Modeling and Design: A Review\.Note:*Archives of Computational Methods in Engineering*, 33:5223–5260doi:10\.1007/s11831\-025\-10448\-9Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- Phan and Luong \(2024\)Q\. M\. Phan and N\. H\. LuongEfficient multi\-fidelity neural architecture search with zero\-cost proxy\-guided local search\.Note:In*Proceedings of the Genetic and Evolutionary Computation Conference \(GECCO ’24\)*, ACMdoi:10\.1145/3638529\.3654027Cited by:[§3\.4](https://arxiv.org/html/2609.23023#S3.SS4.p1.1)\.
- Qinet al\.\(2026\)Z\. Qin, C\. Wang, J\. Guo, H\. Qin, P\. Shen, and W\. JiaBLADE: Brainstorming LLMs as algorithm designer through evolution for online subset selection\.Note:*Knowledge\-Based Systems*, 343:115890doi:10\.1016/j\.knosys\.2026\.115890Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Raissiet al\.\(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. KarniadakisPhysics\-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Note:*Journal of Computational Physics*, 378:686–707Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1)\.
- Rathoreet al\.\(2024\)P\. Rathore, W\. Lei, Z\. Frangella, L\. Lu, and M\. UdellChallenges in training PINNs: A loss landscape perspective\.Note:In*Proceedings of the 41st International Conference on Machine Learning*, PMLR 235:42159–42191Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p4.1)\.
- Renet al\.\(2026\)Z\. Ren, S\. Zhou, D\. Liu, and Q\. LiuPINNAS: Physics\-Informed Neural Network Architecture Search for Solving Integral Equations\.Note:*Journal of Computer Science and Technology*, 41\(2\):593–608doi:10\.1007/s11390\-025\-5433\-9Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1),[§1](https://arxiv.org/html/2609.23023#S1.p4.1)\.
- Romera\-Paredeset al\.\(2024\)B\. Romera\-Paredes, M\. Barekatain, A\. Novikov, M\. Balog, M\. P\. Kumar, E\. Dupont, F\. J\. R\. Ruiz, J\. S\. Ellenberg, P\. Wang, O\. Fawzi, P\. Kohli, and A\. FawziMathematical discoveries from program search with large language models\.Note:*Nature*, 625:468–475doi:10\.1038/s41586\-023\-06924\-6Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p1.1)\.
- Sukumar and Srivastava \(2022\)N\. Sukumar and A\. SrivastavaExact imposition of boundary conditions with distance functions in physics\-informed deep neural networks\.Note:*Computer Methods in Applied Mechanics and Engineering*, 389:114333Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Thawonet al\.\(2026\)I\. Thawon, D\. Vo, T\. Q\. Bui, K\. Rattanamongkhonkun, C\. Chamroon, N\. Tippayawong, Y\. Mona, R\. Wanison, and P\. SuttakulPhysics\-Informed Neural Networks: Current Progress and Challenges in Computational Solid and Structural Mechanics\.Note:*Computer Modeling in Engineering & Sciences*, 146\(2\):2doi:10\.32604/cmes\.2026\.077044Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- Toscanoet al\.\(2025\)J\. D\. Toscano, V\. Oommen, A\. J\. Varghese, Z\. Zou, N\. A\. Daryakenari, C\. Wu, and G\. E\. KarniadakisFrom PINNs to PIKANs: Recent advances in physics\-informed machine learning\.Note:*Machine Learning for Computational Science and Engineering*, 1:15doi:10\.1007/s44379\-025\-00015\-1Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p2.1)\.
- van Stein and Bäck \(2025\)N\. van Stein and T\. BäckLLaMEA: A large language model evolutionary algorithm for automatically generating metaheuristics\.Note:*IEEE Transactions on Evolutionary Computation*, 29\(2\):331–345doi:10\.1109/TEVC\.2024\.3497793Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- van Steinet al\.\(2025\)N\. van Stein, D\. Vermetten, and T\. BäckIn\-the\-loop hyper\-parameter optimization for LLM\-based automated design of heuristics\.Note:*ACM Transactions on Evolutionary Learning and Optimization*doi:10\.1145/3731567Cited by:[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- Wang and Hu \(2025\)R\. Wang and D\. HuNon\-uniform random walk for adaptive sampling\.Note:*Journal of Computational Physics*, 538:114160doi:10\.1016/j\.jcp\.2025\.114160Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Wanget al\.\(2024\)S\. Wang, S\. Sankaran, and P\. PerdikarisRespecting causality is all you need for training physics\-informed neural networks\.Note:*Computer Methods in Applied Mechanics and Engineering*, 421:116813Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p3.1)\.
- Wanget al\.\(2021a\)S\. Wang, Y\. Teng, and P\. PerdikarisUnderstanding and mitigating gradient flow pathologies in physics\-informed neural networks\.Note:*SIAM Journal on Scientific Computing*, 43\(5\):A3055–A3081Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p2.1)\.
- Wanget al\.\(2021b\)S\. Wang, H\. Wang, and P\. PerdikarisOn the eigenvector bias of Fourier feature networks: From regression to solving multi\-scale PDEs with physics\-informed neural networks\.Note:*Computer Methods in Applied Mechanics and Engineering*, 384:113938Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Wanget al\.\(2022\)S\. Wang, X\. Yu, and P\. PerdikarisWhen and why PINNs fail to train: A neural tangent kernel perspective\.Note:*Journal of Computational Physics*, 449:110768Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p2.1)\.
- Wanget al\.\(2023\)Y\. Wang, X\. Han, C\. Chang, D\. Zha, U\. Braga\-Neto, and X\. HuAuto\-PINN: Understanding and optimizing physics\-informed neural architecture\.Note:*NeurIPS 2023 Workshop on AI for Science: From Theory to Practice*Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p4.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p1.1)\.
- Wang and Zhong \(2024\)Y\. Wang and L\. ZhongNAS\-PINN: Neural architecture search\-guided physics\-informed neural network for solving PDEs\.Note:*Journal of Computational Physics*, 496:112603Cited by:[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p2.1)\.
- Wanget al\.\(2026\)Y\. Wang, L\. Zhang, Y\. Guo, and L\. WuTowards automated optimization algorithm design with LLM: An exploratory study in multi\-objective weather routing\.Note:*Transportation Research Part E: Logistics and Transportation Review*, 210:104786doi:10\.1016/j\.tre\.2026\.104786Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Wonet al\.\(2025\)J\. Won, H\. Lee, and J\. LeeA review on multi\-fidelity hyperparameter optimization in machine learning\.Note:*ICT Express*, 11\(2\):245–257doi:10\.1016/j\.icte\.2025\.02\.008Cited by:[§3\.4](https://arxiv.org/html/2609.23023#S3.SS4.p1.1)\.
- Wuet al\.\(2023\)C\. Wu, M\. Zhu, Q\. Tan, Y\. Kartha, and L\. LuA comprehensive study of non\-adaptive and residual\-based adaptive sampling for physics\-informed neural networks\.Note:*Computer Methods in Applied Mechanics and Engineering*, 403:115671Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p3.1)\.
- Wuet al\.\(2024\)H\. Wu, H\. Luo, Y\. Ma, J\. Wang, and M\. LongRoPINN: Region optimized physics\-informed neural networks\.Note:*Advances in Neural Information Processing Systems*, 37Cited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px3.p1.1)\.
- Wuet al\.\(2025\)X\. Wu, S\. Wu, J\. Wu, L\. Feng, and K\. C\. TanEvolutionary computation in the era of large language model: Survey and roadmap\.Note:*IEEE Transactions on Evolutionary Computation*, 29\(2\):534–554doi:10\.1109/TEVC\.2024\.3506731Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Wuwuet al\.\(2025\)Q\. Wuwu, C\. Gao, T\. Chen, Y\. Huang, Y\. Zhang, J\. Wang, J\. Li, H\. Zhou, and S\. ZhangPINNsAgent: Automated PDE Surrogation with Large Language Models\.Note:In*Proceedings of the 42nd International Conference on Machine Learning*, PMLR 267:68143–68165Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p4.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p4.1),[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px3.p1.1)\.
- Yanget al\.\(2026\)X\. Yang, M\. Ran, and H\. PuA novel adaptive gradient\-enhanced physics\-informed neural networks for solving partial differential equations\.Note:*Computers & Mathematics with Applications*, 213:242–256doi:10\.1016/j\.camwa\.2026\.04\.026Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p2.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Yaoet al\.\(2025\)S\. Yao, F\. Liu, X\. Lin, Z\. Lu, Z\. Wang, and Q\. ZhangMulti\-Objective Evolution of Heuristic Using Large Language Model\.Note:*Proceedings of the AAAI Conference on Artificial Intelligence*, 39\(25\):27144–27152doi:10\.1609/aaai\.v39i25\.34922Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- Yeet al\.\(2024\)H\. Ye, J\. Wang, Z\. Cao, F\. Berto, C\. Hua, H\. Kim, J\. Park, and G\. SongReEvo: Large language models as hyper\-heuristics with reflective evolution\.Note:*Advances in Neural Information Processing Systems*, 37Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.3](https://arxiv.org/html/2609.23023#S2.SS3.p2.1)\.
- Yinet al\.\(2025\)H\. Yin, A\. V\. Kononova, T\. Bäck, and N\. van SteinControlling the mutation in large language models for the efficient evolution of algorithms\.Note:In*Applications of Evolutionary Computation \(EvoApplications 2025\)*, pages 403–417, SpringerCited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Yu and Liu \(2025\)H\. Yu and J\. LiuAutomatically optimizing heuristics for robust scale\-free network design via large language models\.Note:*Scientific Reports*, 15:41030doi:10\.1038/s41598\-025\-25031\-2Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Yu and Liu \(2026\)H\. Yu and J\. LiuExploring automated algorithm design synergizing large language models and evolutionary algorithms: Survey and insights\.Note:*Evolutionary Computation*, 34\(2\):167–194doi:10\.1162/EVCO\.a\.370Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Yuet al\.\(2022\)J\. Yu, L\. Lu, X\. Meng, and G\. E\. KarniadakisGradient\-enhanced physics\-informed neural networks for forward and inverse PDE problems\.Note:*Computer Methods in Applied Mechanics and Engineering*, 393:114823Cited by:[§2\.1](https://arxiv.org/html/2609.23023#S2.SS1.p5.1)\.
- Yuet al\.\(2020\)K\. Yu, C\. Sciuto, M\. Jaggi, C\. Musat, and M\. SalzmannEvaluating the search phase of neural architecture search\.Note:*International Conference on Learning Representations*Cited by:[§3\.4](https://arxiv.org/html/2609.23023#S3.SS4.p2.1)\.
- Zhang and Yang \(2024\)B\. Zhang and C\. YangDiscovering physics\-informed neural networks model for solving partial differential equations through evolutionary computation\.Note:*Swarm and Evolutionary Computation*, 88:101589Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1),[§2\.2](https://arxiv.org/html/2609.23023#S2.SS2.p2.1)\.
- Zhang and Ding \(2025\)J\. Zhang and C\. DingSimple yet effective adaptive activation functions for physics\-informed neural networks\.Note:*Computer Physics Communications*, 307:109428doi:10\.1016/j\.cpc\.2024\.109428Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Zhang and Deng \(2025\)Y\. Zhang and L\. DengLoss\-driven dynamic weight and residual transformation in physics\-informed neural network\.Note:*Neural Networks*, 191:107841doi:10\.1016/j\.neunet\.2025\.107841Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p1.1),[§1](https://arxiv.org/html/2609.23023#S1.p3.1)\.
- Zhaoet al\.\(2024\)Z\. Zhao, X\. Ding, and B\. A\. PrakashPINNsFormer: A transformer\-based framework for physics\-informed neural networks\.Note:*International Conference on Learning Representations*Cited by:[§4\.1](https://arxiv.org/html/2609.23023#S4.SS1.SSS0.Px3.p1.1)\.
- Zhonget al\.\(2025\)R\. Zhong, A\. G\. Hussien, J\. Yu, and M\. MunetomoLLMOA: A novel large language model assisted hyper\-heuristic optimization algorithm\.Note:*Advanced Engineering Informatics*, 64:103042doi:10\.1016/j\.aei\.2024\.103042Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p7.1)\.
- Zhouet al\.\(2025\)W\. Zhou, H\. Song, and X\. ChuAutomated design for physics\-informed modeling with convolutional neural networks\.Note:*Communications Physics*, 8:509doi:10\.1038/s42005\-025\-02414\-5Cited by:[§1](https://arxiv.org/html/2609.23023#S1.p4.1)\.

## Appendix ABenchmark PDE Definitions

This appendix summarizes the governing equations, computational domains, physical parameters, initial and boundary conditions, and reference solutions of the 25 PDE problems used in our experiments\. The formulations follow the corresponding benchmark settings used for training and evaluation\.

We use

Δu=∑i=1d∂2u∂xi2,∂nu=∇u⋅𝒏,\\Delta u=\\sum\_\{i=1\}^\{d\}\\frac\{\\partial^\{2\}u\}\{\\partial x\_\{i\}^\{2\}\},\\qquad\\partial\_\{n\}u=\\nabla u\\cdot\\boldsymbol\{n\},
where𝒏\\boldsymbol\{n\}denotes the outward unit normal\. Unless otherwise specified, all parameters take the values given below\.

### A\.1One\-dimensional Burgers equation \(burgers\_1d\)

The governing equation is

ut\+u​ux−ν​ux​x=0,\(x,t\)∈\[−1,1\]×\[0,1\],u\_\{t\}\+uu\_\{x\}\-\\nu u\_\{xx\}=0,\\qquad\(x,t\)\\in\[\-1,1\]\\times\[0,1\],
with viscosity

ν=0\.01π\.\\nu=\\frac\{0\.01\}\{\\pi\}\.
The initial condition is

u⁡\(x,0\)=−sin⁡\(π​x\),u\(x,0\)=\-\\sin\(\\pi x\),
and homogeneous Dirichlet boundary conditions are imposed at both spatial boundaries:

u⁡\(−1,t\)=u⁡\(1,t\)=0\.u\(\-1,t\)=u\(1,t\)=0\.

### A\.2Two\-dimensional Burgers system \(burgers\_2d\)

Let

𝒖=\(u,v\)𝖳\.\\boldsymbol\{u\}=\(u,v\)^\{\\mathsf\{T\}\}\.
The coupled two\-dimensional Burgers system is

ut\+u​ux\+v​uy−ν⁡\(ux​x\+uy​y\)\\displaystyle u\_\{t\}\+uu\_\{x\}\+vu\_\{y\}\-\\nu\(u\_\{xx\}\+u\_\{yy\}\)=0,\\displaystyle=0,vt\+u​vx\+v​vy−ν⁡\(vx​x\+vy​y\)\\displaystyle v\_\{t\}\+uv\_\{x\}\+vv\_\{y\}\-\\nu\(v\_\{xx\}\+v\_\{yy\}\)=0,\\displaystyle=0,
on

\(x,y,t\)∈\[0,4\]2×\[0,1\],\(x,y,t\)\\in\[0,4\]^\{2\}\\times\[0,1\],
with

The initial fields are prescribed by the benchmark data:

u⁡\(x,y,0\)=u0​\(x,y\),v⁡\(x,y,0\)=v0​\(x,y\),u\(x,y,0\)=u\_\{0\}\(x,y\),\\qquad v\(x,y,0\)=v\_\{0\}\(x,y\),
whereu0u\_\{0\}andv0v\_\{0\}are interpolated fromburgers2d\_init\_u\_0\.datandburgers2d\_init\_v\_0\.dat, respectively\.

Both components are periodic in the two spatial directions:

u⁡\(0,y,t\)\\displaystyle u\(0,y,t\)=u⁡\(4,y,t\),\\displaystyle=u\(4,y,t\),v⁡\(0,y,t\)\\displaystyle v\(0,y,t\)=v⁡\(4,y,t\),\\displaystyle=v\(4,y,t\),u⁡\(x,0,t\)\\displaystyle u\(x,0,t\)=u⁡\(x,4,t\),\\displaystyle=u\(x,4,t\),v⁡\(x,0,t\)\\displaystyle v\(x,0,t\)=v⁡\(x,4,t\)\.\\displaystyle=v\(x,4,t\)\.

### A\.3Gray–Scott reaction–diffusion system \(gray\_scott\_2d\)

The governing equations are

ut\\displaystyle u\_\{t\}=ε1​Δ​u\+b⁡\(1−u\)−u​v2,\\displaystyle=\\varepsilon\_\{1\}\\Delta u\+b\(1\-u\)\-uv^\{2\},vt\\displaystyle v\_\{t\}=ε2​Δ​v−d​v\+u​v2,\\displaystyle=\\varepsilon\_\{2\}\\Delta v\-dv\+uv^\{2\},
on

\(x,y,t\)∈\[−1,1\]2×\[0,200\],\(x,y,t\)\\in\[\-1,1\]^\{2\}\\times\[0,200\],
with

b=0\.04,d=0\.1,ε1=10−5,ε2=5×10−6\.b=0\.04,\\qquad d=0\.1,\\qquad\\varepsilon\_\{1\}=10^\{\-5\},\\qquad\\varepsilon\_\{2\}=5\\times 10^\{\-6\}\.
The initial conditions are

u⁡\(x,y,0\)=1−exp⁡\[−80​\(\(x\+0\.05\)2\+\(y\+0\.02\)2\)\],u\(x,y,0\)=1\-\\exp\\left\[\-80\\left\(\(x\+0\.05\)^\{2\}\+\(y\+0\.02\)^\{2\}\\right\)\\right\],
and

v⁡\(x,y,0\)=exp⁡\[−80​\(\(x−0\.05\)2\+\(y−0\.02\)2\)\]\.v\(x,y,0\)=\\exp\\left\[\-80\\left\(\(x\-0\.05\)^\{2\}\+\(y\-0\.02\)^\{2\}\\right\)\\right\]\.
Following the benchmark protocol, no additional spatial boundary loss is imposed\.

### A\.4Heat equation on a perforated domain \(heat\_2d\_complex\_geometry\)

The governing equation is

ut−Δ​u=0,\(x,y,t\)∈Ω×\[0,3\]\.u\_\{t\}\-\\Delta u=0,\\qquad\(x,y,t\)\\in\\Omega\\times\[0,3\]\.
The spatial domain is the rectangle

R=\[−8,8\]×\[−12,12\],R=\[\-8,8\]\\times\[\-12,12\],
with 17 circular holes removed:

Ω=R∖\(⋃𝒄∈𝒞1B1​\(𝒄\)∪⋃𝒄∈𝒞0\.4B0\.4​\(𝒄\)\),\\Omega=R\\setminus\\left\(\\bigcup\_\{\\boldsymbol\{c\}\\in\\mathcal\{C\}\_\{1\}\}B\_\{1\}\(\\boldsymbol\{c\}\)\\cup\\bigcup\_\{\\boldsymbol\{c\}\\in\\mathcal\{C\}\_\{0\.4\}\}B\_\{0\.4\}\(\\boldsymbol\{c\}\)\\right\),
where

𝒞1=\{\\displaystyle\\mathcal\{C\}\_\{1\}=\\\{\(−4,−3\),\(4,−3\),\(−4,3\),\(4,3\),\\displaystyle\(\-4,\-3\),\(4,\-3\),\(\-4,3\),\(4,3\),\(−4,−9\),\(4,−9\),\(−4,9\),\(4,9\),\\displaystyle\(\-4,\-9\),\(4,\-9\),\(\-4,9\),\(4,9\),\(0,0\),\(0,6\),\(0,−6\)\},\\displaystyle\(0,0\),\(0,6\),\(0,\-6\)\\\},
and

𝒞0\.4=\{\(−3\.2,−6\),\(−3\.2,6\),\(3\.2,−6\),\(3\.2,6\),\(−3\.2,0\),\(3\.2,0\)\}\.\\mathcal\{C\}\_\{0\.4\}=\\\{\(\-3\.2,\-6\),\(\-3\.2,6\),\(3\.2,\-6\),\(3\.2,6\),\(\-3\.2,0\),\(3\.2,0\)\\\}\.
The initial condition is

Robin boundary conditions are imposed as

∂nu=\{5−u,on the boundaries of the radius\-1 holes,1−u,on the boundaries of the radius\-0\.4 holes,0\.1−u,on the outer rectangular boundary\.\\partial\_\{n\}u=\\begin\{cases\}5\-u,&\\text\{on the boundaries of the radius\-1 holes\},\\\\ 1\-u,&\\text\{on the boundaries of the radius\-0\.4 holes\},\\\\ 0\.1\-u,&\\text\{on the outer rectangular boundary\}\.\\end\{cases\}

### A\.5Long\-time nonlinear heat equation \(heat\_2d\_long\_time\)

The governing equation is

ut−10−3​\(ux​x\+uy​y\)−5​sin⁡\(u2\)​\(1\+2​sin⁡π​t4\)​sin⁡\(4​π​x\)​sin⁡\(2​π​y\)=0,\\displaystyle u\_\{t\}\-10^\{\-3\}\(u\_\{xx\}\+u\_\{yy\}\)\-5\\sin\(u^\{2\}\)\\left\(1\+2\\sin\\frac\{\\pi t\}\{4\}\\right\)\\sin\(4\\pi x\)\\sin\(2\\pi y\)=0,
for

\(x,y,t\)∈\[0,1\]2×\[0,100\]\.\(x,y,t\)\\in\[0,1\]^\{2\}\\times\[0,100\]\.
The initial condition is

u⁡\(x,y,0\)=sin⁡\(4​π​x\)​sin⁡\(3​π​y\)\.u\(x,y,0\)=\\sin\(4\\pi x\)\\sin\(3\\pi y\)\.
Homogeneous Dirichlet boundary conditions are imposed on the spatial boundary:

u⁡\(x,y,t\)=0,\(x,y\)∈∂\[0,1\]2\.u\(x,y,t\)=0,\\qquad\(x,y\)\\in\\partial\[0,1\]^\{2\}\.

### A\.6Multiscale anisotropic heat equation \(heat\_2d\_multiscale\)

The governing equation is

ut−1\(500​π\)2​ux​x−1π2​uy​y=0,u\_\{t\}\-\\frac\{1\}\{\(500\\pi\)^\{2\}\}u\_\{xx\}\-\\frac\{1\}\{\\pi^\{2\}\}u\_\{yy\}=0,
on

\(x,y,t\)∈\[0,1\]2×\[0,5\]\.\(x,y,t\)\\in\[0,1\]^\{2\}\\times\[0,5\]\.
The initial condition is

u⁡\(x,y,0\)=sin⁡\(20​π​x\)​sin⁡\(π​y\)\.u\(x,y,0\)=\\sin\(20\\pi x\)\\sin\(\\pi y\)\.
Homogeneous Dirichlet boundary conditions are imposed:

u⁡\(x,y,t\)=0,\(x,y\)∈∂\[0,1\]2\.u\(x,y,t\)=0,\\qquad\(x,y\)\\in\\partial\[0,1\]^\{2\}\.
The analytical reference solution is

u⋆​\(x,y,t\)=sin⁡\(20​π​x\)​sin⁡\(π​y\)​exp⁡\[−\(\(20​π\)2\(500​π\)2\+π2π2\)​t\]\.u^\{\\star\}\(x,y,t\)=\\sin\(20\\pi x\)\\sin\(\\pi y\)\\exp\\left\[\-\\left\(\\frac\{\(20\\pi\)^\{2\}\}\{\(500\\pi\)^\{2\}\}\+\\frac\{\\pi^\{2\}\}\{\\pi^\{2\}\}\\right\)t\\right\]\.

### A\.7Varying\-coefficient heat equation \(heat\_2d\_varying\_coefficient\)

The governing equation used in the benchmark is

ut−κ⁡\(x,y\)​Δ​u=f⁡\(x,y,t\),u\_\{t\}\-\\kappa\(x,y\)\\Delta u=f\(x,y,t\),
on

\(x,y,t\)∈\[0,1\]2×\[0,5\],\(x,y,t\)\\in\[0,1\]^\{2\}\\times\[0,5\],
whereκ⁡\(x,y\)\\kappa\(x,y\)is the heterogeneous coefficient field stored inheat\_2d\_coef\_256\.datand evaluated by nearest\-neighbour interpolation\.

The source term is

f⁡\(x,y,t\)=200​sin⁡\(π​x\)​sin⁡\(5​π​y\)​sin⁡\(π​t\)\.f\(x,y,t\)=200\\sin\(\\pi x\)\\sin\(5\\pi y\)\\sin\(\\pi t\)\.
The initial condition is

and homogeneous Dirichlet boundary conditions are imposed:

u\|∂\[0,1\]2=0\.u\|\_\{\\partial\[0,1\]^\{2\}\}=0\.

### A\.8Five\-dimensional heat equation \(heat\_5d\)

Let

Ω=\{𝒙∈ℝ5:‖𝒙‖2≤1\}\.\\Omega=\\left\\\{\\boldsymbol\{x\}\\in\\mathbb\{R\}^\{5\}:\\\|\\boldsymbol\{x\}\\\|\_\{2\}\\leq 1\\right\\\}\.
The governing equation is

ut−15​Δ​u=f⁡\(𝒙,t\),\(𝒙,t\)∈Ω×\[0,1\],u\_\{t\}\-\\frac\{1\}\{5\}\\Delta u=f\(\\boldsymbol\{x\},t\),\\qquad\(\\boldsymbol\{x\},t\)\\in\\Omega\\times\[0,1\],
where

f⁡\(𝒙,t\)=−‖𝒙‖225​exp⁡\(‖𝒙‖222\+t\)\.f\(\\boldsymbol\{x\},t\)=\-\\frac\{\\\|\\boldsymbol\{x\}\\\|\_\{2\}^\{2\}\}\{5\}\\exp\\left\(\\frac\{\\\|\\boldsymbol\{x\}\\\|\_\{2\}^\{2\}\}\{2\}\+t\\right\)\.
Define

g⁡\(𝒙,t\)=exp⁡\(‖𝒙‖222\+t\)\.g\(\\boldsymbol\{x\},t\)=\\exp\\left\(\\frac\{\\\|\\boldsymbol\{x\}\\\|\_\{2\}^\{2\}\}\{2\}\+t\\right\)\.
The initial condition is

u⁡\(𝒙,0\)=g⁡\(𝒙,0\),u\(\\boldsymbol\{x\},0\)=g\(\\boldsymbol\{x\},0\),
and the Neumann boundary condition is

∂nu=g⁡\(𝒙,t\),𝒙∈∂Ω\.\\partial\_\{n\}u=g\(\\boldsymbol\{x\},t\),\\qquad\\boldsymbol\{x\}\\in\\partial\\Omega\.
The reference solution is

u⋆​\(𝒙,t\)=g⁡\(𝒙,t\)\.u^\{\\star\}\(\\boldsymbol\{x\},t\)=g\(\\boldsymbol\{x\},t\)\.

### A\.9Inverse heat equation \(heat\_inverse\_2d\)

The unknowns are the stateu⁡\(x,y,t\)u\(x,y,t\)and the spatially varying diffusion coefficienta⁡\(x,y\)a\(x,y\)\. They satisfy

ut−∇⋅\(a∇u\)=f\(x,y,t\),u\_\{t\}\-\\nabla\\cdot\(a\\nabla u\)=f\(x,y,t\),
on

\(x,y,t\)∈\[−1,1\]2×\[0,1\]\.\(x,y,t\)\\in\[\-1,1\]^\{2\}\\times\[0,1\]\.
The manufactured reference fields are

u⋆​\(x,y,t\)=e−t​sin⁡\(π​x\)​sin⁡\(π​y\),u^\{\\star\}\(x,y,t\)=e^\{\-t\}\\sin\(\\pi x\)\\sin\(\\pi y\),
and

a⋆​\(x,y\)=2\+sin⁡\(π​x\)​sin⁡\(π​y\)\.a^\{\\star\}\(x,y\)=2\+\\sin\(\\pi x\)\\sin\(\\pi y\)\.
The corresponding source term is

f\(x,y,t\)=e−t\[\\displaystyle f\(x,y,t\)=e^\{\-t\}\\Big\[\(4​π2−1\)​sin⁡\(π​x\)​sin⁡\(π​y\)\\displaystyle\(4\\pi^\{2\}\-1\)\\sin\(\\pi x\)\\sin\(\\pi y\)\+π2​\(2​sin2⁡\(π​x\)​sin2⁡\(π​y\)−cos2⁡\(π​x\)​sin2⁡\(π​y\)CLOSE\\displaystyle\+\\pi^\{2\}\\big\(2\\sin^\{2\}\(\\pi x\)\\sin^\{2\}\(\\pi y\)\-\\cos^\{2\}\(\\pi x\)\\sin^\{2\}\(\\pi y\)−sin2\(πx\)cos2\(πy\)\)\]\.\\displaystyle\-\\sin^\{2\}\(\\pi x\)\\cos^\{2\}\(\\pi y\)\\big\)\\Big\]\.
Since

on the spatial boundary, the coefficient boundary condition is

a⁡\(x,y\)=2,\(x,y\)∈∂\[−1,1\]2\.a\(x,y\)=2,\\qquad\(x,y\)\\in\\partial\[\-1,1\]^\{2\}\.
The inverse problem additionally uses 2500 noisy state observations:

uobs=u⋆\+ϵ,ϵ∼𝒩⁡\(0,0\.12\)\.u\_\{\\mathrm\{obs\}\}=u^\{\\star\}\+\\epsilon,\\qquad\\epsilon\\sim\\mathcal\{N\}\(0,0\.1^\{2\}\)\.

### A\.10Kuramoto–Sivashinsky equation \(kuramoto\_sivashinsky\_1d\)

The governing equation is

ut\+α​u​ux\+β​ux​x\+γ​ux​x​x​x=0,u\_\{t\}\+\\alpha uu\_\{x\}\+\\beta u\_\{xx\}\+\\gamma u\_\{xxxx\}=0,
on

\(x,t\)∈\[0,2​π\]×\[0,1\],\(x,t\)\\in\[0,2\\pi\]\\times\[0,1\],
where

α=10016,β=100162,γ=100164\.\\alpha=\\frac\{100\}\{16\},\\qquad\\beta=\\frac\{100\}\{16^\{2\}\},\\qquad\\gamma=\\frac\{100\}\{16^\{4\}\}\.
The initial condition is

u⁡\(x,0\)=cos⁡\(x\)​\(1\+sin⁡x\)\.u\(x,0\)=\\cos\(x\)\(1\+\\sin x\)\.
Following the benchmark protocol, the spatial domain is treated as periodic, while only the initial condition is explicitly included in the constraint loss\.

### A\.11Navier–Stokes flow over a backward\-facing step \(navier\_stokes\_2d\_backstep\)

Let

𝒖=\(u,v\)𝖳\.\\boldsymbol\{u\}=\(u,v\)^\{\\mathsf\{T\}\}\.
The steady incompressible Navier–Stokes equations are

\(𝒖⋅∇\)𝒖\+∇p−νΔ𝒖=𝟎,\(\\boldsymbol\{u\}\\cdot\\nabla\)\\boldsymbol\{u\}\+\\nabla p\-\\nu\\Delta\\boldsymbol\{u\}=\\boldsymbol\{0\},
together with

∇⋅𝒖=0,\\nabla\\cdot\\boldsymbol\{u\}=0,
where

The computational domain is

Ω=\(\[0,4\]×\[0,2\]\)∖\(\[0,2\]×\[1,2\]\)\.\\Omega=\(\[0,4\]\\times\[0,2\]\)\\setminus\(\[0,2\]\\times\[1,2\]\)\.
At the inletx=0x=0,0≤y≤10\\leq y\\leq 1,

u⁡\(0,y\)=4​y​\(1−y\),v⁡\(0,y\)=0\.u\(0,y\)=4y\(1\-y\),\\qquad v\(0,y\)=0\.
At the outletx=4x=4,

The no\-slip condition is imposed on the remaining solid walls:

### A\.12Lid\-driven cavity flow \(navier\_stokes\_2d\_classic\)

This problem corresponds to the benchmark casenavier\_stokes\_2d\_C\.

The steady incompressible Navier–Stokes equations are

\(𝒖⋅∇\)𝒖\+∇p−νΔ𝒖=𝟎,\(\\boldsymbol\{u\}\\cdot\\nabla\)\\boldsymbol\{u\}\+\\nabla p\-\\nu\\Delta\\boldsymbol\{u\}=\\boldsymbol\{0\},
and

∇⋅𝒖=0,\\nabla\\cdot\\boldsymbol\{u\}=0,
on

\(x,y\)∈\[0,1\]2\.\(x,y\)\\in\[0,1\]^\{2\}\.
The viscosity is

corresponding to

On the moving top wall,

u⁡\(x,1\)=4​x​\(1−x\),v⁡\(x,1\)=0\.u\(x,1\)=4x\(1\-x\),\\qquad v\(x,1\)=0\.
On the other three walls,

Pressure is fixed by the gauge condition

### A\.13Long\-time incompressible Navier–Stokes equation \(navier\_stokes\_2d\_long\_time\)

The governing equations are

𝒖t\+\(𝒖⋅∇\)𝒖\+∇p−νΔ𝒖=𝒇\(x,y,t\),\\boldsymbol\{u\}\_\{t\}\+\(\\boldsymbol\{u\}\\cdot\\nabla\)\\boldsymbol\{u\}\+\\nabla p\-\\nu\\Delta\\boldsymbol\{u\}=\\boldsymbol\{f\}\(x,y,t\),
and

∇⋅𝒖=0,\\nabla\\cdot\\boldsymbol\{u\}=0,
on

\(x,y,t\)∈\[0,2\]×\[0,1\]×\[0,5\],\(x,y,t\)\\in\[0,2\]\\times\[0,1\]\\times\[0,5\],
with

The forcing term is

𝒇⁡\(x,y,t\)=\(0−sin⁡\(π​x\)​sin⁡\(π​y\)​sin⁡\(π​t\)\)\.\\boldsymbol\{f\}\(x,y,t\)=\\begin\{pmatrix\}0\\\\ \-\\sin\(\\pi x\)\\sin\(\\pi y\)\\sin\(\\pi t\)\\end\{pmatrix\}\.
At the inletx=0x=0,

u⁡\(0,y,t\)=sin⁡\(π​y\)​\[sin⁡\(π​t\)\+sin⁡\(3​π​t\)\+sin⁡\(5​π​t\)\],u\(0,y,t\)=\\sin\(\\pi y\)\\left\[\\sin\(\\pi t\)\+\\sin\(3\\pi t\)\+\\sin\(5\\pi t\)\\right\],
and

At the outletx=2x=2,

On the remaining spatial boundaries,

The initial condition is

u⁡\(x,y,0\)=v⁡\(x,y,0\)=p⁡\(x,y,0\)=0\.u\(x,y,0\)=v\(x,y,0\)=p\(x,y,0\)=0\.

### A\.14Classical Poisson equation on a perforated domain \(poisson\_2d\_classic\)

The governing equation is

Δ​u=0in​Ω,\\Delta u=0\\qquad\\text\{in \}\\Omega,
where

Ω=\[−12,12\]2∖⋃𝒄∈𝒞B0\.1​\(𝒄\),\\Omega=\\left\[\-\\frac\{1\}\{2\},\\frac\{1\}\{2\}\\right\]^\{2\}\\setminus\\bigcup\_\{\\boldsymbol\{c\}\\in\\mathcal\{C\}\}B\_\{0\.1\}\(\\boldsymbol\{c\}\),
and

𝒞=\{\(0\.3,0\.3\),\(−0\.3,0\.3\),\(0\.3,−0\.3\),\(−0\.3,−0\.3\)\}\.\\mathcal\{C\}=\\\{\(0\.3,0\.3\),\(\-0\.3,0\.3\),\(0\.3,\-0\.3\),\(\-0\.3,\-0\.3\)\\\}\.
The Dirichlet boundary conditions are

u=1on the outer rectangular boundary,u=1\\qquad\\text\{on the outer rectangular boundary\},
and

u=0on the four circular boundaries\.u=0\\qquad\\text\{on the four circular boundaries\}\.

### A\.15Screened Poisson–Boltzmann problem \(poisson\_boltzmann\_2d\)

The benchmark considers the screened Poisson equation

−Δ​u\+k2​u=f⁡\(x,y\),\-\\Delta u\+k^\{2\}u=f\(x,y\),
with

The source term is

f⁡\(x,y\)=A⁡\(μ12\+x2\+μ22\+y2\)​sin⁡\(μ1​π​x\)​sin⁡\(μ2​π​y\),f\(x,y\)=A\\left\(\\mu\_\{1\}^\{2\}\+x^\{2\}\+\\mu\_\{2\}^\{2\}\+y^\{2\}\\right\)\\sin\(\\mu\_\{1\}\\pi x\)\\sin\(\\mu\_\{2\}\\pi y\),
where

A=10,μ1=1,μ2=4\.A=10,\\qquad\\mu\_\{1\}=1,\\qquad\\mu\_\{2\}=4\.
The spatial domain is

Ω=\[−1,1\]2∖⋃j=14Brj​\(𝒄j\),\\Omega=\[\-1,1\]^\{2\}\\setminus\\bigcup\_\{j=1\}^\{4\}B\_\{r\_\{j\}\}\(\\boldsymbol\{c\}\_\{j\}\),
where

\(𝒄1,r1\)\\displaystyle\(\\boldsymbol\{c\}\_\{1\},r\_\{1\}\)=\(\(0\.5,0\.5\),0\.2\),\\displaystyle=\(\(0\.5,0\.5\),0\.2\),\(𝒄2,r2\)\\displaystyle\(\\boldsymbol\{c\}\_\{2\},r\_\{2\}\)=\(\(0\.4,−0\.4\),0\.4\),\\displaystyle=\(\(0\.4,\-0\.4\),0\.4\),\(𝒄3,r3\)\\displaystyle\(\\boldsymbol\{c\}\_\{3\},r\_\{3\}\)=\(\(−0\.2,−0\.7\),0\.1\),\\displaystyle=\(\(\-0\.2,\-0\.7\),0\.1\),\(𝒄4,r4\)\\displaystyle\(\\boldsymbol\{c\}\_\{4\},r\_\{4\}\)=\(\(−0\.6,0\.5\),0\.3\)\.\\displaystyle=\(\(\-0\.6,0\.5\),0\.3\)\.
The Dirichlet conditions are

u=0\.2on the outer boundary,u=0\.2\\qquad\\text\{on the outer boundary\},
and

u=1on the circular boundaries\.u=1\\qquad\\text\{on the circular boundaries\}\.

### A\.16Multiscale multi\-area Poisson equation \(poisson\_2d\_many\_area\)

The governing equation is

a⁡\(x,y\)​Δ​u\+f⁡\(x,y\)=0,\(x,y\)∈\[−10,10\]2\.a\(x,y\)\\Delta u\+f\(x,y\)=0,\\qquad\(x,y\)\\in\[\-10,10\]^\{2\}\.
The domain is partitioned into5×55\\times 5rectangular subdomains\. Within subdomainΩi​j\\Omega\_\{ij\},

a⁡\(x,y\)=ai​j,a\(x,y\)=a\_\{ij\},
whereai​ja\_\{ij\}is read frompoisson\_a\_coef\.dat\.

The local forcing is defined by the tabulated coefficients inpoisson\_f\_coef\.dat:

f⁡\(x,y\)=ci​j,00\+∑r=01∑s=01ci​j,r​s​sin⁡\(π​r​ξhx\)​sin⁡\(π​s​ηhy\),\(x,y\)∈Ωi​j,f\(x,y\)=c\_\{ij,00\}\+\\sum\_\{r=0\}^\{1\}\\sum\_\{s=0\}^\{1\}c\_\{ij,rs\}\\sin\\left\(\\frac\{\\pi r\\,\\xi\}\{h\_\{x\}\}\\right\)\\sin\\left\(\\frac\{\\pi s\\,\\eta\}\{h\_\{y\}\}\\right\),\\qquad\(x,y\)\\in\\Omega\_\{ij\},
where\(ξ,η\)\(\\xi,\\eta\)denote the local coordinates withinΩi​j\\Omega\_\{ij\}\.

The Robin boundary condition is

∂nu=−uon​∂\[−10,10\]2\.\\partial\_\{n\}u=\-u\\qquad\\text\{on \}\\partial\[\-10,10\]^\{2\}\.

### A\.17Three\-dimensional Poisson equation on a complex geometry \(poisson\_3d\_complex\_geometry\)

The governing equation is

−μ⁡\(z\)​Δ​u\+κ⁡\(z\)​u=f⁡\(x,y,z\),\-\\mu\(z\)\\Delta u\+\\kappa\(z\)u=f\(x,y,z\),
in a unit cube with four spherical holes:

Ω=\[0,1\]3∖⋃j=14Brj​\(𝒄j\),\\Omega=\[0,1\]^\{3\}\\setminus\\bigcup\_\{j=1\}^\{4\}B\_\{r\_\{j\}\}\(\\boldsymbol\{c\}\_\{j\}\),
where

\(𝒄1,r1\)\\displaystyle\(\\boldsymbol\{c\}\_\{1\},r\_\{1\}\)=\(\(0\.4,0\.3,0\.6\),0\.2\),\\displaystyle=\(\(0\.4,0\.3,0\.6\),0\.2\),\(𝒄2,r2\)\\displaystyle\(\\boldsymbol\{c\}\_\{2\},r\_\{2\}\)=\(\(0\.6,0\.7,0\.6\),0\.2\),\\displaystyle=\(\(0\.6,0\.7,0\.6\),0\.2\),\(𝒄3,r3\)\\displaystyle\(\\boldsymbol\{c\}\_\{3\},r\_\{3\}\)=\(\(0\.2,0\.8,0\.7\),0\.1\),\\displaystyle=\(\(0\.2,0\.8,0\.7\),0\.1\),\(𝒄4,r4\)\\displaystyle\(\\boldsymbol\{c\}\_\{4\},r\_\{4\}\)=\(\(0\.6,0\.2,0\.3\),0\.1\)\.\\displaystyle=\(\(0\.6,0\.2,0\.3\),0\.1\)\.
The coefficients are

and

κ⁡\(z\)=\{82,z<0\.5,102,z≥0\.5\.\\kappa\(z\)=\\begin\{cases\}8^\{2\},&z<0\.5,\\\\ 10^\{2\},&z\\geq 0\.5\.\\end\{cases\}
Let

r2=x2\+y2\+z2\.r^\{2\}=x^\{2\}\+y^\{2\}\+z^\{2\}\.
The forcing term is

f⁡\(x,y,z\)=\\displaystyle f\(x,y,z\)=\{\}20​exp⁡\[sin⁡\(π​x\)\+sin⁡\(10​π​y\)\+sin⁡\(5​π​z\)\]​r2−1r2\+1\\displaystyle 20\\exp\\left\[\\sin\(\\pi x\)\+\\sin\(10\\pi y\)\+\\sin\(5\\pi z\)\\right\]\\frac\{r^\{2\}\-1\}\{r^\{2\}\+1\}\+100​sin⁡\(π​x\)​sin⁡\(10​π​y\)​sin⁡\(5​π​z\)\.\\displaystyle\+100\\sin\(\\pi x\)\\sin\(10\\pi y\)\\sin\(5\\pi z\)\.
A homogeneous Neumann condition is imposed on the complete boundary:

∂nu=0on​∂Ω\.\\partial\_\{n\}u=0\\qquad\\text\{on \}\\partial\\Omega\.

### A\.18Five\-dimensional Poisson equation \(poisson\_5d\)

The governing equation is

−Δ​u=π24​∑i=15sin⁡\(π​xi2\),𝒙∈\[0,1\]5\.\-\\Delta u=\\frac\{\\pi^\{2\}\}\{4\}\\sum\_\{i=1\}^\{5\}\\sin\\left\(\\frac\{\\pi x\_\{i\}\}\{2\}\\right\),\\qquad\\boldsymbol\{x\}\\in\[0,1\]^\{5\}\.
The analytical reference solution is

u⋆​\(𝒙\)=∑i=15sin⁡\(π​xi2\)\.u^\{\\star\}\(\\boldsymbol\{x\}\)=\\sum\_\{i=1\}^\{5\}\\sin\\left\(\\frac\{\\pi x\_\{i\}\}\{2\}\\right\)\.
The Dirichlet boundary condition is

u=u⋆on​∂\[0,1\]5\.u=u^\{\\star\}\\qquad\\text\{on \}\\partial\[0,1\]^\{5\}\.

### A\.19Inverse Poisson equation \(poisson\_inverse\_2d\)

The unknowns areu⁡\(x,y\)u\(x,y\)and the spatially varying diffusion coefficienta⁡\(x,y\)a\(x,y\), satisfying

−∇⋅\(a∇u\)=f\(x,y\),\(x,y\)∈\[0,1\]2\.\-\\nabla\\cdot\(a\\nabla u\)=f\(x,y\),\\qquad\(x,y\)\\in\[0,1\]^\{2\}\.
The manufactured reference fields are

u⋆​\(x,y\)=sin⁡\(π​x\)​sin⁡\(π​y\),u^\{\\star\}\(x,y\)=\\sin\(\\pi x\)\\sin\(\\pi y\),
and

a⋆​\(x,y\)=11\+x2\+y2\+\(x−1\)2\+\(y−1\)2\.a^\{\\star\}\(x,y\)=\\frac\{1\}\{1\+x^\{2\}\+y^\{2\}\+\(x\-1\)^\{2\}\+\(y\-1\)^\{2\}\}\.
Writing

a=a⋆​\(x,y\),a=a^\{\\star\}\(x,y\),
the corresponding source term is

f⁡\(x,y\)=\\displaystyle f\(x,y\)=\{\}2​π2​a​sin⁡\(π​x\)​sin⁡\(π​y\)\\displaystyle 2\\pi^\{2\}a\\sin\(\\pi x\)\\sin\(\\pi y\)\+2πa2\[\(2x−1\)cos\(πx\)sin\(πy\)\\displaystyle\+2\\pi a^\{2\}\\Big\[\(2x\-1\)\\cos\(\\pi x\)\\sin\(\\pi y\)\+\(2y−1\)sin\(πx\)cos\(πy\)\]\.\\displaystyle\+\(2y\-1\)\\sin\(\\pi x\)\\cos\(\\pi y\)\\Big\]\.
The diffusion coefficient is prescribed on the spatial boundary:

a=a⋆on​∂\[0,1\]2\.a=a^\{\\star\}\\qquad\\text\{on \}\\partial\[0,1\]^\{2\}\.
The inverse problem uses a50×5050\\times 50grid of noisy state observations:

uobs=u⋆\+ϵ,ϵ∼𝒩⁡\(0,0\.12\)\.u\_\{\\mathrm\{obs\}\}=u^\{\\star\}\+\\epsilon,\\qquad\\epsilon\\sim\\mathcal\{N\}\(0,0\.1^\{2\}\)\.

### A\.20One\-dimensional wave equation \(wave\_1d\)

The governing equation is

ut​t−4​ux​x=0,\(x,t\)∈\[0,1\]×\[0,1\]\.u\_\{tt\}\-4u\_\{xx\}=0,\\qquad\(x,t\)\\in\[0,1\]\\times\[0,1\]\.
The initial conditions are

u⁡\(x,0\)=sin⁡\(π​x\)\+12​sin⁡\(4​π​x\),u\(x,0\)=\\sin\(\\pi x\)\+\\frac\{1\}\{2\}\\sin\(4\\pi x\),
and

The spatial boundary condition is

u⁡\(0,t\)=u⁡\(1,t\)=0\.u\(0,t\)=u\(1,t\)=0\.
The analytical reference solution is

u⋆​\(x,t\)=sin⁡\(π​x\)​cos⁡\(2​π​t\)\+12​sin⁡\(4​π​x\)​cos⁡\(8​π​t\)\.u^\{\\star\}\(x,t\)=\\sin\(\\pi x\)\\cos\(2\\pi t\)\+\\frac\{1\}\{2\}\\sin\(4\\pi x\)\\cos\(8\\pi t\)\.

### A\.21Heterogeneous two\-dimensional wave equation \(wave\_2d\_heterogeneous\)

The governing equation is

Δ​u−1c⁡\(x,y\)​ut​t=0,\\Delta u\-\\frac\{1\}\{c\(x,y\)\}u\_\{tt\}=0,
or equivalently,

ut​t=c⁡\(x,y\)​Δ​u,u\_\{tt\}=c\(x,y\)\\Delta u,
on

\(x,y,t\)∈\[−1,1\]2×\[0,5\]\.\(x,y,t\)\\in\[\-1,1\]^\{2\}\\times\[0,5\]\.
The heterogeneous coefficientc⁡\(x,y\)c\(x,y\)is defined by interpolation of the Gaussian\-random\-field realization stored indarcy\_2d\_coef\_256\.dat\.

The initial conditions are

u⁡\(x,y,0\)=exp⁡\[−\(x\+0\.5\)2\+y22​\(0\.3\)2\],u\(x,y,0\)=\\exp\\left\[\-\\frac\{\(x\+0\.5\)^\{2\}\+y^\{2\}\}\{2\(0\.3\)^\{2\}\}\\right\],
and

ut​\(x,y,0\)=0\.u\_\{t\}\(x,y,0\)=0\.
A homogeneous Neumann condition is imposed on the spatial boundary:

∂nu=0on​∂\[−1,1\]2\.\\partial\_\{n\}u=0\\qquad\\text\{on \}\\partial\[\-1,1\]^\{2\}\.

### A\.22One\-dimensional Allen–Cahn equation \(allen\_cahn\_1d\)

The governing equation is

ut−ε2​ux​x\+u3−u=0,u\_\{t\}\-\\varepsilon^\{2\}u\_\{xx\}\+u^\{3\}\-u=0,
on

\(x,t\)∈\[−1,1\]×\[0,0\.25\],\(x,t\)\\in\[\-1,1\]\\times\[0,0\.25\],
with

The initial condition is

u⁡\(x,0\)=tanh⁡\(x2​ε\)\.u\(x,0\)=\\tanh\\left\(\\frac\{x\}\{\\sqrt\{2\}\\varepsilon\}\\right\)\.
The Dirichlet boundary condition is obtained from the stationary reference solution:

u⁡\(±1,t\)=tanh⁡\(±12​ε\)\.u\(\\pm 1,t\)=\\tanh\\left\(\\frac\{\\pm 1\}\{\\sqrt\{2\}\\varepsilon\}\\right\)\.
The analytical reference solution is

u⋆​\(x,t\)=tanh⁡\(x2​ε\)\.u^\{\\star\}\(x,t\)=\\tanh\\left\(\\frac\{x\}\{\\sqrt\{2\}\\varepsilon\}\\right\)\.

### A\.23Variable\-coefficient Darcy equation \(darcy\_flow\_2d\)

The governing equation is

−∇⋅\(a\(x,y\)∇u\)=f\(x,y\),\(x,y\)∈\[0,1\]2\.\-\\nabla\\cdot\\left\(a\(x,y\)\\nabla u\\right\)=f\(x,y\),\\qquad\(x,y\)\\in\[0,1\]^\{2\}\.
The heterogeneous coefficient is

a⁡\(x,y\)=2\+12​sin⁡\(2​π​x\)​sin⁡\(3​π​y\)\.a\(x,y\)=2\+\\frac\{1\}\{2\}\\sin\(2\\pi x\)\\sin\(3\\pi y\)\.
The manufactured reference solution is

u⋆​\(x,y\)=sin⁡\(π​x\)​sin⁡\(π​y\)\.u^\{\\star\}\(x,y\)=\\sin\(\\pi x\)\\sin\(\\pi y\)\.
The source term is defined consistently as

f\(x,y\)=−∇⋅\(a\(x,y\)∇u⋆\(x,y\)\)\.f\(x,y\)=\-\\nabla\\cdot\\left\(a\(x,y\)\\nabla u^\{\\star\}\(x,y\)\\right\)\.
Equivalently,

f=−\(ax​ux⋆\+ay​uy⋆\+a​Δ​u⋆\)\.f=\-\\left\(a\_\{x\}u\_\{x\}^\{\\star\}\+a\_\{y\}u\_\{y\}^\{\\star\}\+a\\Delta u^\{\\star\}\\right\)\.
The homogeneous Dirichlet boundary condition is

u=0on​∂\[0,1\]2\.u=0\\qquad\\text\{on \}\\partial\[0,1\]^\{2\}\.

### A\.24Two\-dimensional shallow\-water system \(shallow\_water\_2d\)

Lethhdenote the water depth, and letqxq\_\{x\}andqyq\_\{y\}denote the two momentum components\.

The conservative shallow\-water equations are

ht\+\(qx\)x\+\(qy\)y=0,h\_\{t\}\+\(q\_\{x\}\)\_\{x\}\+\(q\_\{y\}\)\_\{y\}=0,
\(qx\)t\+\(qx2h\+g​h22\)x\+\(qx​qyh\)y=0,\(q\_\{x\}\)\_\{t\}\+\\left\(\\frac\{q\_\{x\}^\{2\}\}\{h\}\+\\frac\{gh^\{2\}\}\{2\}\\right\)\_\{x\}\+\\left\(\\frac\{q\_\{x\}q\_\{y\}\}\{h\}\\right\)\_\{y\}=0,
and

\(qy\)t\+\(qx​qyh\)x\+\(qy2h\+g​h22\)y=0,\(q\_\{y\}\)\_\{t\}\+\\left\(\\frac\{q\_\{x\}q\_\{y\}\}\{h\}\\right\)\_\{x\}\+\\left\(\\frac\{q\_\{y\}^\{2\}\}\{h\}\+\\frac\{gh^\{2\}\}\{2\}\\right\)\_\{y\}=0,
on

\(x,y,t\)∈\[0,1\]2×\[0,0\.1\],\(x,y,t\)\\in\[0,1\]^\{2\}\\times\[0,0\.1\],
with gravitational acceleration

The initial state is spatially uniform:

\(h,qx,qy\)​\(x,y,0\)=\(1,0\.1,−0\.05\)\.\(h,q\_\{x\},q\_\{y\}\)\(x,y,0\)=\(1,0\.1,\-0\.05\)\.
All three variables are periodic in both spatial directions:

𝑼⁡\(0,y,t\)=𝑼⁡\(1,y,t\),\\boldsymbol\{U\}\(0,y,t\)=\\boldsymbol\{U\}\(1,y,t\),
and

𝑼⁡\(x,0,t\)=𝑼⁡\(x,1,t\),\\boldsymbol\{U\}\(x,0,t\)=\\boldsymbol\{U\}\(x,1,t\),
where

𝑼=\(h,qx,qy\)𝖳\.\\boldsymbol\{U\}=\(h,q\_\{x\},q\_\{y\}\)^\{\\mathsf\{T\}\}\.
Water\-depth positivity is enforced through

h≥hmin,hmin=10−3\.h\\geq h\_\{\\min\},\\qquad h\_\{\\min\}=10^\{\-3\}\.
The reference solution is the constant state

𝑼⋆=\(1,0\.1,−0\.05\)𝖳\.\\boldsymbol\{U\}^\{\\star\}=\(1,0\.1,\-0\.05\)^\{\\mathsf\{T\}\}\.

### A\.25Kovasznay flow \(kovasznay\_flow\_2d\)

The steady incompressible Navier–Stokes equations are

\(𝒖⋅∇\)𝒖\+∇p−νΔ𝒖=𝟎,\(\\boldsymbol\{u\}\\cdot\\nabla\)\\boldsymbol\{u\}\+\\nabla p\-\\nu\\Delta\\boldsymbol\{u\}=\\boldsymbol\{0\},
together with

∇⋅𝒖=0\.\\nabla\\cdot\\boldsymbol\{u\}=0\.
The computational domain is

\(x,y\)∈\[−0\.5,1\]×\[−0\.5,1\.5\]\.\(x,y\)\\in\[\-0\.5,1\]\\times\[\-0\.5,1\.5\]\.
The Reynolds number and viscosity are

Re=40,ν=140\.\\mathrm\{Re\}=40,\\qquad\\nu=\\frac\{1\}\{40\}\.
Define

λ=Re2−Re24\+4​π2\.\\lambda=\\frac\{\\mathrm\{Re\}\}\{2\}\-\\sqrt\{\\frac\{\\mathrm\{Re\}^\{2\}\}\{4\}\+4\\pi^\{2\}\}\.
The analytical reference solution is

u⋆​\(x,y\)=1−eλ​x​cos⁡\(2​π​y\),u^\{\\star\}\(x,y\)=1\-e^\{\\lambda x\}\\cos\(2\\pi y\),
v⋆​\(x,y\)=λ2​π​eλ​x​sin⁡\(2​π​y\),v^\{\\star\}\(x,y\)=\\frac\{\\lambda\}\{2\\pi\}e^\{\\lambda x\}\\sin\(2\\pi y\),
and

p⋆​\(x,y\)=12​\(1−e2​λ​x\)\.p^\{\\star\}\(x,y\)=\\frac\{1\}\{2\}\\left\(1\-e^\{2\\lambda x\}\\right\)\.
The analytical values are imposed on the rectangular boundary:

\(u,v,p\)=\(u⋆,v⋆,p⋆\)on​∂Ω\.\(u,v,p\)=\(u^\{\\star\},v^\{\\star\},p^\{\\star\}\)\\qquad\\text\{on \}\\partial\\Omega\.
Pressure is additionally anchored by

### A\.26Note on the excluded PINNacle Wave2d\-MS problem

Among the PINNacle problems, we exclude the 2D multi\-scale long\-time wave equation \(Wave2d\-MS\) because its published mathematical specification is internally inconsistent\. PINNacle defines the governing equation as

ut​t−\(ux​x\+a2​uy​y\)=0,a=2,u\_\{tt\}\-\(u\_\{xx\}\+a^\{2\}u\_\{yy\}\)=0,\\qquad a=\\sqrt\{2\},
while the reported reference solution is

sinh⁡\(π​x\)​sinh⁡\(π​y\)​cos⁡\(3​π​t\)\.\\sinh\(\\pi x\)\\sinh\(\\pi y\)\\cos\(\\sqrt\{3\}\\pi t\)\.
Direct substitution gives

−6​π2​u⋆≠0\.\-6\\pi^\{2\}u^\{\\star\}\\neq 0\.
showing that the reported reference solution does not satisfy the stated governing equation\.

We further observed that the reference solution in the released PINNacle implementation differs from the one reported in the paper and, under the released default parameters, is also inconsistent with the stated PDE\. To avoid evaluating PINN methods against an ambiguous PDE–reference pair, we therefore exclude Wave2d\-MS from our benchmark rather than modifying its original definition\. All other PINNacle problems used in this work follow the corresponding benchmark formulations described above\.

## Appendix BImplementation Details and Mechanistic Evidence for LLM Candidate Generation

### B\.1LLM Inputs, Outputs, and Information Boundaries

The LLM in PINNsForge does not directly generate arbitrary executable code\. Instead, each call proposes only one complete structuredAlgorithmSpec\. A controller is responsible for subsequent parsing, validation, instantiation, and real PINN training, so the LLM cannot bypass the predefined search space or directly modify the experimental protocol\.

The information provided to the LLM at different search stages is summarized below\.

ContextGeneration 0Generationg≥1g\\geq 1PurposeTarget PDE description✓\\checkmark✓\\checkmarkProvides the governing equation, computational domain, IC/BC, and problem characteristicsPrior knowledgeK𝒫K\_\{\\mathcal\{P\}\}✓\\checkmark✓\\checkmarkProvides PINN design rationale relevant to the current PDEAllowed search spaceΩ\\Omega✓\\checkmark✓\\checkmarkRestricts the network, sampling, constraint, loss, and optimization components that may be generatedTraining / resource budget✓\\checkmark✓\\checkmarkRestricts parameter count, sampling scale, and training budgetRole / search\-strategy condition✓\\checkmark✓\\checkmarkSpecifies the design focus of the current candidateGlobal Top\-3AlgorithmSpec–✓\\checkmarkProvides validated high\-performing parent configurationsExecution evidence of parents–✓\\checkmarkProvides training evidence including PDE residuals, IC/BC errors, optimization stability, and convergence behaviorRun\-scoped posterior summary–✓\\checkmarkSummarizes successes, failures, and design experience observed within the current independent search runAt Generation 0, no execution history exists for the current PDE, so initial candidates are independently generated using only PDE information, relevant priors, the allowed search space, and initialization roles\. Starting from Generation 1, the Global Top\-3 parent configurations and their real execution evidence are added to the context to support the Refine, Physics\-Guided, Architecture\-Guided, Synthesis, and Novelty strategies defined in the main text\. Sampling temperatures for the different strategies are provided in Appendix C\.2; the role definitions are not repeated here\.

The LLM is not allowed to modify the target PDE, computational domain, initial conditions, boundary conditions, evaluation protocol, or experimental budget\. Its output can describe only the PINN algorithm components permitted withinΩ\\Omega, including network representation, network architecture, activation functions, constraint handling, sampling strategy, loss handling, optimization strategy, and training control\. Candidates are normalized before deduplication and training; normalizedAlgorithmSpecs that are exactly identical are treated as the same configuration, avoiding duplicate candidates caused solely by field ordering, expansion of default values, or differences in ineffective representations\.

### B\.2Information Isolation in Ablation Experiments

Section 4\.3 of the main text defines three major ablated variants\. To ensure that differences among these experiments arise only from the designated information channels, we keep the number of candidates, training budget, search space, and LLM generation protocol unchanged while controlling the context accessible to the LLM\.

Table B\.1: Information accessible to the LLM under different ablation settings\.

VariantPrior KnowledgeRun\-scoped PosteriorExecution Feedback to LLMParent\-based Evolutionw/o Knowledge Base×\\times×\\times✓\\checkmark✓\\checkmarkw/o Execution Feedback✓\\checkmark×\\times×\\times✓\\checkmarkw/o Evolutionary Search✓\\checkmark✓\\checkmark✓\\checkmark×\\timesPINNsForge✓\\checkmark✓\\checkmark✓\\checkmark✓\\checkmarkInw/o Knowledge Base, the LLM no longer accesses the retrieved PINN priors or the run\-specific relevant knowledge maintained through this knowledge channel\.

Inw/o Execution Feedback, historical candidate structures remain visible, but theiryg,ky\_\{g,k\}, PDE residuals, constraint errors, optimization diagnostics, rankings, and other training outcomes are hidden\. No posterior depending on these execution results is formed\. The LLM therefore cannot use candidates’ actual execution behavior to guide subsequent modifications\.

Inw/o Evolutionary Search, the LLM can still use aggregated historical execution experience, but it does not receive specific historicalAlgorithmSpecs or Global Top\-3 parent configurations\. Subsequent candidates are therefore generated independently and do not inherit, modify, or recombine the concrete configurations of existing high\-performing candidates\. This setting corresponds to independent LLM search without parent inheritance under the same candidate\-evaluation budget\.

### B\.3Deterministic Validation, Normalization, and Candidate Repair

The LLM is responsible only for proposing candidates; whether a candidate is admitted to training is determined entirely by a deterministic controller\. Every returnedAlgorithmSpecpasses through

Parsing→Normalization→Validation→Instantiation→Training\.\\text\{Parsing\}\\rightarrow\\text\{Normalization\}\\rightarrow\\text\{Validation\}\\rightarrow\\text\{Instantiation\}\\rightarrow\\text\{Training\}\.
The validation stage checks at least the following conditions:

- •whether the structure of theAlgorithmSpecis complete;
- •whether all components belong to the allowed search spaceΩ\\Omega;
- •whether different algorithm components are mutually compatible;
- •whether the number of network parameters satisfies the budget;
- •whether the sampling budget satisfies the constraints;
- •whether the training iteration count strictly matches the current low\- or high\-fidelity training budget\.

Invalid candidates do not enter training directly\. Instead, they undergo structured repair within a fixed repair budget\. Repair may address only missing fields, illegal values, or component incompatibilities; it cannot modify the target PDE𝒫\\mathcal\{P\}, enlarge the search spaceΩ\\Omega, or relax the experimental budget\. If the required number of valid candidates still cannot be obtained within the allowed generation and repair process, low\-fidelity search terminates according to the stopping rule in Appendix C\.4\.

Candidates that pass validation are reconstructed and trained according to the experimental protocol, and their actual execution results form\(yg,k,Eg,k\)\(y\_\{g,k\},E\_\{g,k\}\)for subsequent search\. Natural\-language or structured suggestions produced by the LLM are therefore not themselves considered valid candidates; onlyAlgorithmSpecs that pass deterministic validation and complete real PINN training are admitted to the history and the Global Top\-3\.

### B\.4Execution Feedback and Cross\-Generation Inheritance Evidence on Burgers\-1D

The main text already uses search trajectories to show that multi\-generation search continues to discover better candidates after Generation 0; complete trajectories for all 25 PDEs are provided in Appendix D\. Here, rather than repeating the Burgers equation definition, initialization roles, or generation\-wise global\-best values, we retain only two concrete pieces of evidence illustratinghow real training feedback changes subsequent designandhow the final configuration emerges through multi\-generation parent inheritance\.

Below, we useGg\-Ckas a readable label for candidateag,ka\_\{g,k\}; for example,G6\-C4corresponds toa6,4a\_\{6,4\}\.

##### B\.4\.1 Representative Parent—Child Modification

For parentG6\-C4, real training shows that the network and hard constraints remain stable, but a clear PDE\-residual hotspot persists near the shock\. The subsequent candidate therefore preserves the validated stable network and constraint mechanisms, concentrates its main modification on local adaptive sampling, and simultaneously explores an alternative optimization configuration\.

Table B\.2: Representative execution\-feedback\-driven modification on Burgers\-1D\.

ComponentParent ConfigurationExecution EvidenceDesign DecisionChild ConfigurationArchitectureModified MLP,6×1286\\times 128; no residual connectionsCompleted the full 1,000\-step training without NaN or solution collapsePreserve\.The current network architecture is already stableKeep Modified MLP,6×1286\\times 128ConstraintExact hard transformburgers\_1d\_initial\_dirichlet\_exactIC error=2\.25×10−8=2\.25\\times 10^\{\-8\}, BC error=8\.74×10−8=8\.74\\times 10^\{\-8\}Preserve\.Do not reallocate search emphasis to constraints that are already well satisfiedKeep the exact hard transformSampling10,000 Sobol interior points \+ RAD; add 500 each time; cap at 15,000; residual exponent 2\.5Persistent residual hotspot atx∈\[−0\.24,0\.24\]x\\in\[\-0\.24,0\.24\],t∈\[0\.15,0\.39\]t\\in\[0\.15,0\.39\]Modify\.Strengthen local refinement while preserving global Sobol coverageAdd 600 each time; raise the cap to 16,000; keep the remaining RAD settings unchangedOptimizationAdam 800 steps, lr1\.0×10−31\.0\\times 10^\{\-3\}, one\-cycle; L\-BFGS 200 steps, history 100Optimization remains numerically stable, but clear residual error remains at the end of trainingExplore\.Treat optimization changes as accompanying exploration rather than misclassifying a stable optimizer as a failureIncrease Adam lr to1\.2×10−31\.2\\times 10^\{\-3\}, switch to step decay; L\-BFGS history 120The actual training results of the parent and child are:

MetricParent G6\-C4Child G7\-C4Solution MSE5\.71×10−45\.71\\times 10^\{\-4\}2\.21×10−42\.21\\times 10^\{\-4\}PDE residual6\.05×10−26\.05\\times 10^\{\-2\}2\.14×10−22\.14\\times 10^\{\-2\}IC error2\.25×10−82\.25\\times 10^\{\-8\}1\.25×10−81\.25\\times 10^\{\-8\}BC error8\.74×10−88\.74\\times 10^\{\-8\}2\.53×10−82\.53\\times 10^\{\-8\}This example is not intended to attribute the performance improvement to any single independent component\. The parent and child use independent random initializations, and both the sampling and optimization configurations change simultaneously\. Instead, the example illustrates how execution feedback helps the system distinguish between “mechanisms that are already stable and should be preserved” and “components that still require modification,” and use this distinction to form the next\-generation candidate\.

##### B\.4\.2 Cross\-Generation Evolutionary Lineage of the Final Candidate

PINNsForge records the parent source of every evolutionary candidate\. Single\-parent refinement can be written as

ag,p→ag\+1,c,a\_\{g,p\}\\rightarrow a\_\{g\+1,c\},
Multi\-parent synthesis instead inherits and recombines design mechanisms from multiple Global Top\-3 candidates\. Here, refinement and synthesis are not predefined field\-wise mutation or crossover operations; which designs are preserved, modified, or combined is determined by the LLM from the parentAlgorithmSpecs, the target PDE, and execution evidence\.

Table B\.3: Cross\-generation evolutionary lineage of the final Burgers\-1D candidate\.

GenerationCandidateParent\(s\)Evolution ModeMain Inherited / Modified Components0G0\-C2–Independent Initialization4×\\times64 modified MLP,tanh, soft constraints, Sobol \+ RAD, Adam \+ L\-BFGS1G1\-C3G0\-C2, G0\-C6, G0\-C4Multi\-Parent SynthesisCombine RAD sampling, a deeper network topology, and optimization control; expand to a 6×\\times96 modified MLP2G2\-C1G1\-C3Single\-Parent RefinementPreserve the 6×\\times96 network; increase the RAD exponent and interior\-point budget3G3\-C4G2\-C1, G2\-C4Multi\-Parent SynthesisIntroduce an exact hard IC/BC transform; use fixed weighting; increase the RAD cap4G4\-C4G3\-C4, G3\-C2Multi\-Parent SynthesisPreserve hard constraints and Sobol \+ RAD; widen the network to 6×\\times128; add one\-cycle Adam and L\-BFGS recovery5G5\-C4G4\-C4, G4\-C3Multi\-Parent SynthesisAdjust the Adam schedule, RAD capacity, and L\-BFGS history6G6\-C4G5\-C4, G5\-C1, G5\-C3Multi\-Parent SynthesisPreserve hard constraints,tanh, Sobol \+ RAD, and Adam→\\rightarrowL\-BFGS; use a stable 6×\\times128 modified MLP7G7\-C4G6\-C4Single\-Parent RefinementStrengthen RAD refinement; adjust the Adam schedule and L\-BFGS history8G8\-C4G7\-C1, G7\-C2, G7\-C4Multi\-Parent SynthesisCombine network capacity, sampling, and training strategies from multiple G7 parents; use a 5×\\times192 modified MLP and retune RAD and L\-BFGS settingsFinal HFG8\-C4Low\-fidelity source: G8\-C4Independent High\-Fidelity TrainingKeep the G8\-C4AlgorithmSpecunchanged and independently train for 10,000 optimizer steps from a new random initializationThis lineage shows that the final configuration is not produced directly by a single independent LLM call, but instead continually inherits, modifies, and synthesizes design mechanisms that have been validated through real PINN training over multiple generations\. Together with the search trajectories in the main text and Appendix D, this example provides a traceable mechanistic chain fromexecution evidence→\\rightarrowtargeted modification→\\rightarrowparent inheritance→\\rightarrowfinal configuration\.

## Appendix CExperimental Details

This appendix provides additional experimental details for reproducing the PINNsForge results\. The notation for candidates, execution results, search history, and Global Top\-3 follows Section 3 and Appendix B\.

### C\.1Search Configuration

PINNsForge starts with 8 candidates in Generation 0 and generates 5 candidates in each subsequent generation\. Generation indices start from \(g=0\)\. The search completes at least Generations 0–3 and at most Generations 0–9, corresponding to a minimum of 4 and a maximum of 10 completed generations\. After each completed generation, the controller maintains the distinct all\-history Global Top\-3 defined in Section 3\.3\.3\.

SettingValueGeneration\-0 candidates8Candidates per later generation5Min\. / max\. completed generations4 / 10Parent populationGlobal Top\-3Low\-fidelity training1K stepsHigh\-fidelity training10K stepsMaximum trainable parameters2MMaximum training points65,536All normal and escape generations count toward the 10\-generation maximum\.

### C\.2LLM Sampling

Candidate generation uses a fixed role\-conditioned temperature schedule\. Lower temperatures are assigned to conservative search directions, while higher temperatures encourage exploratory designs\.

Generation\-0 directionTemperatureRobust reference0\.20Constraint / causality0\.35Localized sampling0\.40Representation0\.45Residual / loss0\.45Optimization0\.40Capacity / topology0\.50Cross\-component novelty0\.60Evolution strategyNormalEscapeRefine0\.200\.35Physics\-Guided0\.360\.51Architecture\-Guided0\.420\.57Synthesis0\.470\.62Novelty0\.580\.73When stagnation triggers an escape generation, the temperature of each evolution strategy is increased by 0\.15, capped at 0\.85\. The Reflection Agent uses a temperature of 0\.1\.top\_pfollows the provider default, and no provider\-side sampling seed is supplied\.

### C\.3Computational Environment

Experiments are conducted on a Linux server using NVIDIA GeForce RTX 4090 GPUs\. Each PDE search runs on a single GPU using FP32 arithmetic\.

ComponentConfigurationGPUNVIDIA GeForce RTX 4090NVIDIA Driver550\.163\.01CUDA12\.4Python3\.12PyTorch2\.6\.0\+cu124PrecisionFP32Additional software dependencies and exact environment specifications are provided with the implementation\.

### C\.4Deterministic Stagnation Detection, Escape Generation, and Termination

All search\-control decisions in this section are computed from the low\-fidelity candidate comparison metric defined in Section 3\.3\.2\. For candidateag,ka\_\{g,k\}, let \(y\_\{g,k\}\) denote its low\-fidelity comparison value\. In the reported experiments, \(y\_\{g,k\}\) is the solution MSE computed exclusively on the fixed search\-validation reference set\. The disjoint final\-test reference set is never accessed by search\-control decisions, candidate ranking, or termination\. Lower values are better\. For each completed generation \(g\), define the best valid candidate value in that generation as

bg=mink⁡yg,k,b\_\{g\}=\\min\_\{k\}y\_\{g,k\},
where the minimum is taken over valid candidates in Generationgg\. Define the all\-history best value through Generationggas

Bg=min0≤t≤g⁡bt\.B\_\{g\}=\\min\_\{0\\leq t\\leq g\}b\_\{t\}\.
When Generationggcontains at least three valid candidates, define the median of its three best candidate values as

mg=Median⁡\(Top3⁡\{yg,k\}k\)\.m\_\{g\}=\\operatorname\{Median\}\\left\(\\operatorname\{Top3\}\\\{y\_\{g,k\}\\\}\_\{k\}\\right\)\.
The median comparison value of the distinct all\-history Global Top\-3 archive after completing Generationggis defined as

Ag=Median⁡\(\{y⁡\(a\):a∈GlobalDistinctTop3⁡\(ℋg\)\}\)\.A\_\{g\}=\\operatorname\{Median\}\\left\(\\left\\\{y\(a\):a\\in\\operatorname\{GlobalDistinctTop3\}\(\\mathcal\{H\}\_\{g\}\)\\right\\\}\\right\)\.
This quantity is used once at least three distinct valid candidates exist in the archive; before that point, archive\-median criteria are disabled\.

The relative improvement from an earlier valuexxto a later valueyyis defined as

I⁡\(x,y\)=\{max⁡\(0,x−y\|x\|\),x≠0​and both​x,y​are finite,0,otherwise\.I\(x,y\)=\\begin\{cases\}\\max\\left\(0,\\dfrac\{x\-y\}\{\|x\|\}\\right\),&x\\neq 0\\text\{ and both \}x,y\\text\{ are finite\},\\\\\[8\.0pt\] 0,&\\text\{otherwise\}\.\\end\{cases\}
With this definition, largerI⁡\(x,y\)I\(x,y\)indicates a larger reduction in the comparison metric\.

Stagnation detection begins after Generation 3 has completed, i\.e\., after four completed generations \(Generations 0–3\)\. For a normal generationg≥3g\\geq 3, an escape generation is scheduled when at least one of the following conditions holds\.

1. 1\.Global\-best stagnation\.The three most recent generation\-to\-generation improvements of the all\-history best are all below 1%, and the cumulative improvement over the same window is below 2%:

I⁡\(Bg−3,Bg−2\)<0\.01,I⁡\(Bg−2,Bg−1\)<0\.01,I\(B\_\{g\-3\},B\_\{g\-2\}\)<0\.01,\\qquad I\(B\_\{g\-2\},B\_\{g\-1\}\)<0\.01,
I⁡\(Bg−1,Bg\)<0\.01,I⁡\(Bg−3,Bg\)<0\.02\.I\(B\_\{g\-1\},B\_\{g\}\)<0\.01,\\qquad I\(B\_\{g\-3\},B\_\{g\}\)<0\.02\.
1. 2\.Elite\-archive stagnation\.Across the last three completed generations, no previously unseen normalizedAlgorithmSpecenters the distinct all\-history Global Top\-3 archive\.
2. 3\.Population stagnation\.When each generation in the relevant window contains at least three valid candidates, the median of the three best candidates improves by less than 2% in two consecutive transitions:

I⁡\(mg−2,mg−1\)<0\.02,I⁡\(mg−1,mg\)<0\.02\.I\(m\_\{g\-2\},m\_\{g\-1\}\)<0\.02,\\qquad I\(m\_\{g\-1\},m\_\{g\}\)<0\.02\.
If stagnation is detected after Generationgstallg\_\{\\mathrm\{stall\}\}, and the 10\-generation limit has not yet been reached, the next generation is designated as the escape generation:

gesc=gstall\+1\.g\_\{\\mathrm\{esc\}\}=g\_\{\\mathrm\{stall\}\}\+1\.
The escape generation uses the same five semantic evolution strategies and the same 1,000\-step low\-fidelity training budget as a normal evolution generation, but increases exploration through the following deterministic controls:

- •each evolution\-strategy temperature is increased by 0\.15, capped at 0\.85;
- •every generated candidate must differ from its parent context in at least two major executable modules;
- •exact reproductions of parent configurations and non\-substantive numerical perturbations are rejected;
- •previously failed mechanisms are not repeated unless the candidate explicitly addresses the observed failure reason;
- •Synthesis and Novelty are prompted to use greater parent or mechanism diversity\.

The normal evolution\-strategy temperature vector is

𝝉normal=\(0\.20,0\.36,0\.42,0\.47,0\.58\),\\boldsymbol\{\\tau\}\_\{\\mathrm\{normal\}\}=\(0\.20,\\,0\.36,\\,0\.42,\\,0\.47,\\,0\.58\),
and the corresponding escape\-generation vector is

𝝉escape=\(0\.35,0\.51,0\.57,0\.62,0\.73\)\.\\boldsymbol\{\\tau\}\_\{\\mathrm\{escape\}\}=\(0\.35,\\,0\.51,\\,0\.57,\\,0\.62,\\,0\.73\)\.
After completing the escape generationgescg\_\{\\mathrm\{esc\}\}, the escape is considered successful if at least one of the following conditions holds\.

1. 1\.New best candidate\.The best candidate in the escape generation improves on the pre\-escape all\-history best by at least 2%:

I⁡\(Bgstall,bgesc\)≥0\.02\.I\\left\(B\_\{g\_\{\\mathrm\{stall\}\}\},b\_\{g\_\{\\mathrm\{esc\}\}\}\\right\)\\geq 0\.02\.
1. 2\.New elite configuration\.At least one previously unseen normalizedAlgorithmSpecfrom Generationgescg\_\{\\mathrm\{esc\}\}enters the distinct all\-history Global Top\-3 archive\.
2. 3\.Elite\-archive improvement\.The median comparison value of the distinct Global Top\-3 improves by at least 2%:

I⁡\(Agstall,Agesc\)≥0\.02\.I\\left\(A\_\{g\_\{\\mathrm\{stall\}\}\},A\_\{g\_\{\\mathrm\{esc\}\}\}\\right\)\\geq 0\.02\.
If the escape succeeds, the stagnation history accumulated beforegescg\_\{\\mathrm\{esc\}\}is cleared and normal evolution resumes\. If none of the three escape\-success conditions is satisfied, the low\-fidelity search terminates immediately and the current Global Top\-3 is passed to the high\-fidelity evaluation stage\.

The low\-fidelity search also terminates when either of the following conditions holds:

- •the search has completed the maximum of 10 generations, i\.e\., Generation 9 has finished; or
- •candidate generation and the bounded repair process cannot produce the required number of valid candidates for the next generation\.

Normal and escape generations both count toward the 10\-generation maximum\. Because generation indices run from 0 to 9 and stagnation is first assessed after Generation 3, the earliest possible escape generation is Generation 4\.

## Appendix DComplete Search Trajectories

Table D\.1 reports the complete single\-run search trajectories for all 25 PDE problems\. All MSE values used to describe the search trajectory are computed on the search\-validation reference set\. Generation\-0 Best denotes the lowest validation MSE among the initial candidates, Best LF denotes the lowest validation MSE found during the entire low\-fidelity search, and Final HF denotes the lowest high\-fidelity validation MSE among the independently trained Global Top\-3 candidates\. The final\-test reference set is not used in constructing these search trajectories\. The improvement ratios are computed as

G0→Best LF=MSEG0MSEBest​LF,\\text\{G0\}\\rightarrow\\text\{Best LF\}=\\frac\{\\mathrm\{MSE\}\_\{\\mathrm\{G0\}\}\}\{\\mathrm\{MSE\}\_\{\\mathrm\{Best\\ LF\}\}\},
G0→Final HF=MSEG0MSEFinal​HF\.\\text\{G0\}\\rightarrow\\text\{Final HF\}=\\frac\{\\mathrm\{MSE\}\_\{\\mathrm\{G0\}\}\}\{\\mathrm\{MSE\}\_\{\\mathrm\{Final\\ HF\}\}\}\.
Table D\.1: Complete single\-run search trajectories across 25 PDE problems\.

PDEGeneration\-0 Best MSEBest LF MSEFinal HF MSEG0→\\rightarrowBest LFG0→\\rightarrowFinal HFburgers\_1d3\.64E\-032\.05E\-04 \(G8\)6\.51E\-0517\.74×\\times55\.82×\\timesburgers\_2d2\.49E\-012\.13E\-01 \(G5\)1\.78E\-011\.17×\\times1\.39×\\timesgray\_scott\_2d4\.37E\-034\.03E\-03 \(G9\)4\.07E\-031\.08×\\times1\.07×\\timesheat\_2d\_complex\_geometry1\.28E\-023\.26E\-03 \(G8\)4\.63E\-043\.92×\\times27\.64×\\timesheat\_2d\_long\_time1\.13E\+001\.10E\+00 \(G5\)1\.08E\+001\.03×\\times1\.05×\\timesheat\_2d\_multiscale3\.01E\-029\.140E\-05 \(G9\)1\.014E\-05330\.08×\\times2,975\.09×\\timesheat\_2d\_varying\_coefficient1\.12E\-022\.85E\-03 \(G9\)2\.18E\-033\.93×\\times5\.13×\\timesheat\_5d1\.31E\-053\.85E\-07 \(G9\)8\.99E\-0933\.90×\\times1,452\.78×\\timesheat\_inverse\_2d1\.36E\-025\.73E\-03 \(G6\)3\.68E\-032\.36×\\times3\.68×\\timeskuramoto\_sivashinsky\_1d1\.34E\+001\.05E\+00 \(G7\)1\.16E\+001\.28×\\times1\.16×\\timesnavier\_stokes\_2d\_backstep2\.64E\-037\.68E\-04 \(G9\)2\.30E\-043\.44×\\times11\.48×\\timesnavier\_stokes\_2d\_classic1\.86E\-031\.61E\-05 \(G9\)2\.35E\-06115\.60×\\times792\.09×\\timesnavier\_stokes\_2d\_long\_time5\.09E\+025\.08E\+02 \(G7\)4\.97E\+021\.00×\\times1\.02×\\timespoisson\_2d\_classic9\.77E\-028\.94E\-02 \(G8\)1\.64E\-021\.09×\\times5\.95×\\timespoisson\_boltzmann\_2d1\.04E\-011\.23E\-02 \(G9\)6\.03E\-038\.45×\\times17\.16×\\timespoisson\_2d\_many\_area3\.27E\+002\.79E\+00 \(G7\)1\.97E\+001\.17×\\times1\.66×\\timespoisson\_3d\_complex\_geometry8\.96E\-034\.48E\-03 \(G8\)8\.83E\-042\.00×\\times10\.14×\\timespoisson\_5d6\.44E\-061\.01E\-07 \(G9\)5\.65E\-0863\.58×\\times114\.03×\\timespoisson\_inverse\_2d6\.22E\-034\.32E\-04 \(G8\)1\.45E\-0414\.39×\\times42\.80×\\timeswave\_1d7\.43E\-025\.57E\-02 \(G6\)3\.15E\-021\.34×\\times2\.36×\\timeswave\_2d\_heterogeneous7\.85E\-023\.85E\-02 \(G7\)2\.96E\-022\.03×\\times2\.65×\\timesallen\_cahn\_1d6\.37E\-081\.08E\-08 \(G9\)1\.20E\-095\.89×\\times53\.1×\\timesdarcy\_flow\_2d6\.94E\-093\.14E\-10 \(G8\)1\.53E\-1122\.11×\\times453\.59×\\timesshallow\_water\_2d2\.73E\-092\.33E\-10 \(G9\)9\.83E\-1111\.68×\\times27\.71×\\timeskovasznay\_flow\_2d2\.00E\-054\.95E\-08 \(G9\)2\.59E\-09404\.28×\\times7,722\.18×\\times