DGNet: Discrete Green Networks for Data-Efficient Learning of Spatiotemporal PDEs
Authors: Yingjie Tan
Quanming Yao
Yaqing Wang
Venue: ICLR 2026 Poster (2026)
Keywords: Partial Differential Equations
Data-efficient learning
Graph Neural Networks
Physics-Informed Machine Learning
Generalization
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: A data-efficient method for learning spatiotemporal PDEs that generalizes to unseen source terms.

Abstract: Spatiotemporal partial differential equations (PDEs) underpin a wide range of scientific and engineering applications. Neural PDE solvers offer a promising alternative to classical numerical methods. However, existing approaches typically require large numbers of training trajectories, while high-fidelity PDE data are expensive to generate. Under limited data, their performance degrades substantially, highlighting their low data efficiency. A key reason is that PDE dynamics embody strong structural inductive biases that are not explicitly encoded in neural architectures, forcing models to learn fundamental physical structure from data. A particularly salient manifestation of this inefficiency is poor generalization to unseen source terms. In this work, we revisit Green’s function theory—a cornerstone of PDE theory—as a principled source of structural inductive bias for PDE learning. Based on this insight, we propose DGNet, a discrete Green network for data-efficient learning of spatiotemporal PDEs. The key idea is to transform the Green’s function into a graph-based discrete formulation, and embed the superposition principle into the hybrid physics–neural architecture which reduces the burden of learning physical priors from data, thereby improving sample efficiency. Across diverse spatiotemporal PDE scenarios, DGNet consistently achieves state-of-the-art accuracy using only tens of training trajectories. Moreover, it exhibits robust zero-shot generalization to unseen source terms, serving as a stress test that highlights its data-efficient structural design.

Published as a conference paper at ICLR 2026
DGNET: DISCRETE GREEN NETWORKS FOR DATA- EFFICIENT LEARNING OF SPATIOTEMPORAL PDEs
Yingjie Tan
Quanming Yao
Qiuzhen College
Department of Electronic Engineering Tsinghua University qyaoaaâ singhua edu.ct cn
Tsinghua University
tany)23 @ma is. t.s.nnguanaaaa edu. cn
Yaqing Wang"
Beijing Institute Mathematical Sciences and Applications wangya ng@bi msaa.
ABSTRACT
Spatiotemporal partial differential equations (PDEs) underpin wide range of sci- entific and engincering applications. Neural PDE solvers offer promising al- ternative to classical numerical methods. However, cxisting approaches typically require large numbers of training trajectories, while high-fidelity PDE data are expensive to generate. Under limited data, their performance degrades substan- tially, highlighting their low data efficiency. key reason that PDE dynamics embody strong structural inductive biases that are not explicitly encoded in neural architectures, forcing models learn fundamental physical structure from data. A particularly salient manifestation of this inefficiency is poor generalization to unseen source terms. In this work, we revisit Green's function theory-a corner- stone of PDE theory-as a principled source of structural inductive bias for PDE learning. Based on this insight, we propose DGNet, discrete Green network for data-efficient learning of spatiotemporal PDEs. The key idea is to transform the Green's function into a graph-based discrete formulation, and embed the super- position principle into the hybrid physics nnysral architecture, which reduces the burden of learning physical priors from data, thereby improving sample efficiency. Across diverse spatiotemporal PDE scenarios, DGNet consistently achieves state- of-the-art accuracy using only tens of training trajectories. Moreover, it exhibits robust zero-shot generalization unseen source terms, serving stress test that highlights its data-efficient structural design.
INTRODUCTION
Spatiotemporal partial differential equations (PDEs) form the foundation of modeling dynamical systems across science and engineering, governing phenomena in fluid dynamics (Hirsch) 2007). weather forecasting (Lynch] 2008), molecular dynamics (Lelievre & Stoltz 2016), and energy sys- tems (Reos-Mercado & Borraz-Sánchez 2015). Accurately solving such PDEs is crucial for scien- tific discovery and engineering design. Classical numerical solvers (Anderson) 2002 Evans al.) 2012) can provide reliable solutions but often become computationally prohibitive for large-scale irregular domains. Neural PDE solvers have therefore emerged a promising alternative, aim- ing to approximate solution operators directly from data (Raissi et al. 2019 Brandstetter et al. 2022 Zeng et al. 2025). practice, however, high- high-fidelity PDE data whether collected from real systems or generated through simulation-art are expensive to obtain. Existing neural PDE solvers typically require large numbers of training trajectories to achieve strong performance across the so- lution domain. Under limited data regimes, their accuracy degrades substantially, indicating low data efficiency. As a result, improving the data efficiency of neural PDE learning remains a critical yet underexplored challenge.
"Correspondence author

---

Published as a conference paper at ICLR 2026
A fundamental reason for this limitation lies in the structural nature of PDE dynamics. PDE systems embody strong physical and mathematical inductive biases such as locality, conservation laws, and superposition- that govern how states evolve space and time. Yet most neural architectures do not explicitly encode these structural priors (Li et al. 2021 Lu et al.
2021 Sanchez-Gionza
et al.
2020b), only incorporate them weakly Sanchez-Gionzalee
al. 2020a Horie & Mitsume [2022
Bishnoi et al. 2024). Consequently, models must learn the underlying operator structure and phys- ical constraints directly from data, leading to high sample complexity and reduced data efficiency. This issue becomes particularly pronounced in the presence of source terms f(x, t),, which repre- sent external forcing applied over space and time. Examples include time-varying heat sources in conduction (Hahn & Ozisik 20122), body forces in fluid dynamics (Pope), 2001), ), time-dependent cur- rents in electromagneties (Taflove et al. 2005). Training trajectories typically cover only subset of possible source patterns, existing solvers often struggle to extrapolate to unseen source patterns. Poor generalization to unseen source terms thus serves as a representative manifestation of insuffi- cient inductive bias and data inefficiency in current neural PDE approaches.
To improve data efficiency, it is therefore essential to introduce stronger structural priors that mean- ingfully reduce the learning space and exploration burden of neural solvers. Motivated by this perspective, we revisit Green's 's function theory, a foundational principle in PDE analysis. Green's s function provides principled foundation, it effectively extracts the solution from the direct action of the operator fxx recasting the PDE solution into two explicit components- the homogeneous evolution and the forced response. This decomposition encodes the superposition structure intrinsic to many PDE systems, which provides a principled form of inductive bias that reduces the need to rediscover fundamental behavior from data. Moreover, this formulation naturally separates system evolution from source responses, offering a principled foundation for improved generalization under varying source conditions.
Building on this insight, we develop DGNet for spatiotemporal PDE learning on irregular domains commonly encountered in scientific and engineering applications. The core idea to construct a discrete Green's 'ss formulation update rule. By embedding this discrete operator into a hybrid physics-neura architecture, we combine principled physical priors with learnable corrections from graph neural networks (GNNs). This design reduces the burden of learning fundamental dynamics from data while retaining flexi- bility on complex geometries. Empirically, DGNet demonstrates strong data efficiency and robust generalization to unseen source terms.
graphs that preserves the superposition structure in computable
Our contributions can be summarized as follows:
We attribute the limited data efficiency of neural PDE solvers to insufficient structural inductive bias. We provide the insight that Green's ' function theory offers principled prior, which can substantially reduce the need rediscover fundamental physical structure from data.
We introduce DGNet, discrete Green network. By explicitly embedding the superposition prin- ciple into the hybrid physics-nent neural architecture, DGNet encodes physical structure directly in the model rather than requiring it to be rediscovered from data, thereby promoting data-efficient learning by design.
Empirically, DGNet achieves aftte-of-thh-att accuracy across diverse spatiotemporal PDE systems using only tens of training trajectories, highlighting its data-efficien structural design. Further- more, exhibits robust zero-shot generalization to unseen source terms.
RELATED WORK
Physics-informed neural networks (PINNs). PINNs (Raissi et al., 2019) incorporate PDE residu- als and boundary conditions as soft constraints in the training loss, enabling models to fit solutions at sparsely sampled collocation points. However, PINNs essentially Iearn specific solution under a fixed setup, and thus cannot generalize to domains, parameters, or forcing terms beyond the train- ing distribution. Variants (Yu et al., 2022 Sukumar & Srivastava 2022 Costabal al. 2024) have attempted to improve stability or computational efficiency, yet data requirements grow rapidly when accurate solutions are required across the solution domain.
Neural operators. Neural operator methods directly learn mappings between infinite-dimensoon dimensional function spaces from data. Representative works include DeepONet (Lu 2021) and the Fourier
2

---

Published as a conference paper at ICLR 2026
Neural Operator (FNO) al. 2021 along with many extensions in the spectral or kernel do- mains (uuibas et al. 2021 Li et al. 2023; George et al. 2024 Tran et al. 2023). These models achieve generalization across resolutions and physical parameters. However, they typically incor- porate limited physical structure into the model, and learning operator mappings in function spaces requires substantial amounts of data.
Graph-based PDE solvers. For spatiotemporal PDEs on irregular domains, graphs provide a nat- ural backbone. The Graph Network Simulator (GNS) (Sanchez-Gonzalez et al. (020b) introduced
the encooder-procesor-deeodee framework for physical simulation. Subsequent works enhanced physical consistency by embedding conservation laws (Sanchez-
et al., 2020a Crammer
et al. |020] [Bishnoi et al., (023). incorporating operator blocks (Seo et al. 2019 Horie & Mitsumee
2022
Zeng et al., 2025), or enforcing invariants such as momentum conservation (Bishnoi et al,
2024 BENO (Wang et al. 2024) further drew inspiration from Green's 's functions, but focused
on time-independent elliptic PDEs. Despite the incorporation of various inductive biases, existing graph-based solvers lack a principled structural foundation, still rely on substantial data to learn complex physical dynamics.
Generalization challenges in PDE solvers. Generalization is central to neural PDE solvers. Neural operators such as FNO (Li et al. 2021) generalize across resolutions, DeepONet (Lu et al., 2021) across parameters, and BroGNet (Bishnoi al. 2024) across system sizes. Yet critical gap
mains: none of these approaches can generalize to unseen source terms, which frequently arise in scientific and engineering systems. Such varying sources ubiquitous in practice, including time- varying heat sources in conduction (Hahn & Ózisik, 2012), body forces in fluid dynamics (Pope) 2001), time-dependent currents in electromagnetics (Tatlove al. (2005). and seismic excitations in elasticity (Aki & Richards, 2002).
3 DGNET: DISCRETE GREEN NETWORKS
In this section we introduce DGNet, our proposed model for data-efficient spatiotemporal PDE learning enabled by structural inductive bias. We begin in Section [.I by formalizing the problem setup. In Section 3.2 we introduce Green's function and the superposition principle, which naturally decompose the PDE solution into state evolution and source response. the discrete Green formulation by discretizing both space and time, leading preserves the superposition structure. In Section we describe our physics- neural hybrid operator that combines numerical discretization with GNN-based correction construct high-fidelity solver. Finally, Section 3.3 describes the prediction and training procedure.
Section 3.3] we derive
update rule that
3.1 PROBLEM FORMULATION
We consider physical systems governed by spatiotemporal partial differential equations (PDEs) as
Ou(x, t) Øt
f(x,t), (x, t) E Ω x (0,7),
(1)
where u(x,,) t) E IR denotes the physical state, Cx is a spatial differential operator characterizing the system dynamics, and f(x, t) is the source term.
The evolution uniquely determined by an initial condition (IC) and boundary conditions (BCs). The IC specifies u(x, 0) uo(x). On the boundary 02, primarily consider Dirichlet BCs u(x,1) 1) = (xx, 1) for x E aΩD and Neumann BCs (u(x.t) = ln(x, t) for x E aΩN, with M2p U 0ΩN 02 and dΩp n N 0. Other BC types (c.g., Robin, periodic) are accommodated in the same framework.
on
In many real-world systems, the source f(x, t) varies substantially across space and time and thus plays a decisive role in the solution behavior. In practice, models are typically trained on only a lim- ited subset of possible source patterns, and their performance often degrades when confronted with unseen sources. Generalization to varying source terms serves as a stringent test of data efficiency and structural inductive bias in neural PDE solvers.
3

---

Published as a conference paper at ICLR 2026
3.2 GREEN'S FUNCTION REPRESENTATION
Directly solving the PDE I is often challenging- The Green's function method from PDE the- ory (Arfken al. 2013) provides principled way to characterize the system response. The Green's
function G(x, 1;x', is uniquely determined by
0 -x)
- fx)
G(x, = 6(x x') 0(t - T),
(2)
where 8() is the Dirac delta function. Physically, G(x, 1;, 7) describes the influence observed (x, () due to a unit-strength point source applied at (x', T).
Crucially, the Green's function acts mathematical device to "extract" the solution from the direct action of the operator Lx. Instead of learning Lx [u] an inseparable whole, the Green representation rewrites the PDE solution in terms of propagation kernel defined by Cx and convolution with the source term. According to the superposition principle (Boyce et al. 2017), the complete solution u(x, t) can therefore be expressed as
u(x, 1) 1 G 1; (x(xx) dx' + [(
(,,,,xX,,, (x) dx'd.
Evolution of initial state
Respouse to source term
This representation makes two key points explicit. First, the solution naturally decomposes into the evolution the initial condition and the accumulated response to the source term. Second, it enables unified treatment of diverse and time-varying forcing conditions. In the next section, we describe how this continuous formulation is discretized on graphs, leading to a discrete Green's 's function that forms the foundation of our learning architecture and provides strong structural inductive bias for data-efficient learning.
3.3 DISCRETIZATION
While Green's function offers an elegant continuous representation, it is seldom available in closed form for complex domains, and direct evaluation of the associated integrals is computationally pro- hibitive. We therefore derive discrete Green formulation on graphs that preserves the superposition structure within tractable update rule. In this discrete setting, the Green's 's function naturally acts as propagation operator, evolving the system state while incorporating source effects at each time step.
Graph-based ODE system. Graph-baser based discretization is particularly suitable for spatiotemporal PDEs on irregular meshes, as it provides flexible representation of spatial operators and naturally supports message-passing based corrections. Therefore, in DGNet, we discretize the spatial domain Ω into N nodes { (x....., XxN....XNN-1 and represent it as graph G (V, E), where V is the set of nodes and E is the edge set constructed once using Delaunay triangulation. Each node i e V corresponds to spatial location x,.
At discrete time t,, the physical state values {u(x;, ( form state vector u (u(xo, ) ..... u(xN-1,*))
RN,, and similarly the source term values [f(x,, [*)] W-1 form vector p* (f(x, ()), f(xN-1, (²)) R When the time index k is critical, simply denote them by u and f. The continuous PDE T can then be written in discrete form as
du di
Lu +
where L NxN is the matrix discretization of the spatial operator Lx. For node i, the action of on the state vector u is [Lu]; = Enn [Lks [u]j, where [:lk, denotes the (i.j)-th entry of matrix and Fl denotes the j-th component a vector.
Time integration. We apply midpoint (Crank-Nicolson) discretization M /x+1]: let At 4-11 1g. the time increment is (u++1 n²)/At, and the right-hand side is approximated by a trapezoidal average, yielding
2 (Lu + Lu + [(p + p*+1).
(5)
At
4

---

Published as a conference paper at ICLR 2026
Continuous
Discrete
Mesh information
Problem Fomuuation
diseretization
Graph-based ODE system
physiic
Eepural
Green's Function Representation
Discrete Green's 's function
Operator
System evolution
system state a
Source response
predicted state a
source term f
Figure I: Overview of DGNet architecture. The model centers on a bybrid operator Lohysies + Lueurat, where Lppysics encodes gradient and Laplacian discretizations and Lmeural is GNN-based correction for mesh-induced errors. This operator is integrated into the discrete Green's 's function update (Eq. 7]. which naturally combines system evolution with source-term response and provides structural inductive bias for data-efficient learning.
Rearranging terms yields sparse linear system:
II A L u+1 (I
A L) t At (f* + +
(6)
Discrete Green's function. From Eq. the discrete Green's function is identified
G(Al) (I- 2/4
Thus the single-step update becomes
h+1 G(Ar) (1 2 u² G(At) As (P* p*+1).
This discrete formulation mirrors the superposition principle 33: the next state u++ results from the propagation of the current state under G(At) together with the accumulated response to the source term within the interval At. A full derivation of this update scheme is provided in Appendix [
By embedding structural priors directly into the model, DGNet narrows the space of admissible solutions consistent with physical laws. This reduces the burden data recover solution operator and improves learning efficiency, particularly in limited-trajectory regimes.
Since naively computing the inverse (I - L) ()) is infeasible for large systems. Instead. solve the equivalent sparse linear system, where is already sparse due to local mesh interactions. Importantly, the coefficient matrix (I L depends only the static mesh geometry and thus remains fixed throughout rollout. We therefore adopt a "factorize once, solve many times" strategy: cached and reused
single sparse LU factorization is performed before rollout, and its factors for all time steps via efficient forward/backw:. substitution. This reduces per-step cost to nearly linear in the number of nonzeros, making the discrete Green solver practical for large meshes.
3.4 PHYSICS--NURAI HYBRID OPERATOR
The core of our architecture lies in constructing a high-fidelity operator L for the discrete ODE system. Instead of relying solely on data-driven models or purely hand-crafted discretizations, we propose hybrid approach that combines the best of both worlds:
L Lphysics Lacural-
(8)
Here, Lphysics is built directly from mesh geometry using numerical discretization techniques, ensur- ing consistency with physical laws. Meanwhile, Lneural is a learnable correction term, parameterized by GNN, that compensates for discretization irregular meshes.

---

Published as a conference paper at ICLR 2026
3.4.1 PHYSICS PRIOR OPERATOR LpuYSIcs
Among the various spatial operators appearing in PDEs, the gradient and Laplacian are the most fundamental and widely used components (e.g., in diffusion, prrecttonndiifusion and Poisson- type equations). Accordingly, we construct Lphyiis directly from mesh geometry using estab- lished discretization schemes, namely the Green-Gauss theorem (Löhner) 2008) and the dis- crete Laplace-Beltramt Beltrami operator (Meyer et al. 2003). These methods are widely used in com- putational physies and provide physically consistent priors irregular meshes, ensuring that Lphysics encodes reliable physics knowledge. Accordingly, Lphysics is instantiated as follows:
Gradient operator. phyysiesly myyha 2/02,
if E (i))
midee and, pppyssicslii -- ->ENN()) 2 other elements
ay
are zero. Where [2,] is the control volume, and ma is the scalar projection of unit normal nij, li is the face length.
B;
Laplacian operator. [Lphsiissss
Miss if E N(1)
m
and physsicss [i. ΣxEN()) Wake other elements
zero. With cotangent weights Waj (cot ozj cot B.y) and Voronoi area [s2]|-
Figure 2: Geometric variables for the
discrete operator.
Where Ω2 is the control region of node i, (i) is the first-order neighbor set of node i, au and Big are the two angles opposite edge (i,j). The geometric variables are illustrated in Figure
3.4.2 NEURAL CORRECTION OPERATOR INEURAL
Although Lphysics provides a principled foundation, it fundamentally approximates a continuous op- erator discretized mesh, where truncation errors inherent and inevitable. Even though these errors may diminish on high-resolution or regular meshes, the gap between discrete and continuous dynamics persists. To bridge this gap and recover operator consistency without resorting com- putationally expensive fine meshing, we introduce correction matrix Leeurtl learned by GNN following the Encode-Process Decode paradigm (Sanchez SSanclez-Gonzalez al. 2020b).
Encoder: Node features consist of the spatial coordinates and node type C₁ (interior or boundary), which are mapped to embeddings hς NodeEncoder([x, Ct]). Edge features consist of the
relative displacement Xij X, and distance dip ||x; x; 1/112, mapped to embeddings as C.y- EdgeEncoder(( X.j. d.
message-passing passing neural networks (MPNNs). At Processor. The processor consists of M layers of message
m-th layer, we obtain
Update h," E Message(h," hn" ...)) + h
h m++1
Decoder. After M layers, for each edge (i, j), we concatenate final node embeddings and pass them through a decoder to predicts edge-level corrections:
Leeural = Decoder M h,
(10)
For simplicity, we implement all neural components, including the node encoder, edge encoder, message function, update function, and decoder, multilayer perceptrons (MLPs). Theoretical analysis of the consistency and stability properties of the learned correction operator is provided in Appendix E
3.5 PREDICTION AND TRAINING
With the hybrid operator L Lphysies Lneurals the system state is updated by the discrete Green solver in Eq. 7 yielding predictions k at each global time step 4t-

---

Published as a conference paper at ICLR 2026
To further improve accuracy. we include lightweight residual GNN correction module in the pre- diction path. This design follows prior physics neural hybrid solvers (Meng & Karniadakis, 2020; Wu et al., 2024b Long al. 2023 Zeng et al. 2025), which demonstrate that small residual mod- ules can effectively approximate dynamics not explicitly captured by the physical operator, while preserving interpretability of the operator-based formulation.
For training, the full trajectory (u² (u"....,,"""...... _u-1) is segmented into shorter sub-sequences of length Q < T: (u*0 uu", ,..., -* where Sy denotes the global time index of the q-th step in the sub-sequence (corresponding to time point ls.). The model rolls out from u to predict (ū").. .u..... ) Following prior work (Brandstetter et al., [022] Pfaff et al. 2021] Sanchez- Gonzalez et al., 2020b). we apply the pushforward trick and inject small noise into u during training to alleviate error accumulation and improve robustness. The loss is computed only on the first and last predictions of each sub-sequence:
= - - u²I|" + || u²
The learnable parameters are optimized with Adam using a decayed learning rate schedule.
4 EXPERIMENTS
In this section, we empirically evaluate DGNet across diverse scenarios to assess its data efficiency and generalization capability. We consider multiple spatiotemporal PDE systems and examine per- formance under a limited-trajectory training regime. Section 4.1 describes the experimental setup, including datasets, baselines, and evaluation metrics. Section 4.2 presents quantitative and quali- tative results across different categories of PDE systems. Section 4.3 provides ablation studies to analyze the contribution of each architectural component. The source code is publicly available
4.1 EXPERIMENTAL SETUP
Table 1: Governing equations and physical contexts for experimental scenarios.
Name
Physical Scenario
Mathematical Formulation
Meaning of Source Term
Allen-Cahn
Phase separation Population dynamics
duu 2T [u. 12)
Driving force for separation Logistic growth
Fisher-KPP
0u Vu pu(1 u) 8uu = D. (
Naguuo
Excitable systems
D. + ov)
Excitation- recovvrry coupling
Contaminant Trans-
Channel flow with obsta- cles
8ctu Ve= DV -
Reaction (generration/lleaay)
port
MTTT = RV*T - M(T Tame) + Q(211)
Laser Heat
Laser heating
Moving laser heat source
Datasets. We evaluate DGNet on three categories of spatiotemporal PDE systems: (I) classi- cal equations Cal,, Maaher-KP,, FitzHugh-Nagumo) to validate accuracy on canonical diffusion-, advection-, and reaction-driven dynamics, (2) complex geometric domains (contaminant transport with cylinder, sediment, and irregular obstacles) to test robustness on irregular meshes and boundary conditions, and (3) generalization unseen source terms (laser heat treatment) to evaluate adaptability to novel forcing conditions. Table I summarizes the governing equations, physical con- texts, and the meaning of the source term in each scenario. Detailed dataset parameters (domains, mesh sizes, time steps, boundary conditions, and trajectory splits) are deferred to Appendix B
Baselines. We compare DGNet with representative neural PDE solvers: DeepONet (Lu al. 202 MGN (Pfaff et al., 202T). MP-PDE (Brandstetter et al. 2022). PhyMPGN (Zeng et al. 2005)) and BENO (Wang al. 2024). These cover operator-learnin graph-based, and hybrid paradigms. All baselines are tuned to have comparable numbers of parameters and training costs. Implementation details are provided in Appendix
Evaluation Metrics. We report Mean Squared Error (MSE) and Relative & Error (RNE). RNE is defined as || - 1lll2//|uul/ where á and u denote predicted and ground-truth states. We report log MSE for readability, while the relative ranking remains consistent with raw MSE.
/DGGett

---

Published as a conference paper at ICLR 2026
Notably, all experiments are conducted in limited-date regime, with only tens of training tra- jectories per scenario. contrast, neural PDE solvers are often trained with substantially larger datasets ranging from hundreds to thousands of trajectories- achieve stable performance. We evaluate under this low-data setting to highlight the data efficiency of DGNet.
4.2 PERFORMANCE COMPARISON
Table summarizes results across all experimental scenarios, covering the three categories of spa- tiotemporal PDE systems considered in this paper. Even under the limited-trajectory training regime, DGNet achieves stable and state-of-the-ar performance across all tasks and metrics, highlighting the data efficiency of its structural design. Further analysis of efficiency and scalability is provided in Appendix D and additional experimental results on broader scenarios are presented in Appendix
Table 2: Results across scenarios from three categories of spatiotemporal PDE systems under the limited-trajectrry training. For both MSE (in log scale) and RNE, lower values indicate better performance. The best results are highlighted in bold.
Scenario
Metric
DeepONet
MGN
MP-PDE
BENO
PhyMPGN
DGNet
Allen-
MSE
2.60e-01
2.70e-01
8.52e-01
2.52c+00
5.16e-01
8.75e-03
Cahn
RNE
0.6686
0.6813
1.2109
2.0813
0.9420
0.0188
Fisher-
MSE
3.05e-02
3.66e-03
9.90e-02
6.26e-02
1.50e-02
2.59e-04
KPP
RNE
0.4181
0.1448
0.7530
0.5989
0.9270
0.0238
FitzHuugh
MSE
2.49e-06
3.75e-05
6.464e-06
2.14e-04
1.69e-03
1.18e-07
Nagumo
RNE
0.9745
3,4815
1,4454
8.3106
23,5696
0.0952
MSE
4.44e-02
6.38e-03
9.31e-02
6.76e-02
4.13e-01
1.00c-04
Cylinder
RNE
0.5976
0.7154
0.8644
0.7364
1.8201
0.0196
MSE
3.61e-02
5.94e-03
7.10e-03
1.07e-01
2.00e-01
4.60c-04
Sediments
RNE
0.4759
0.6103
0.6673
0.8180
1.1186
0.0282
Complex
MSE
5.33e-02
7.79e-03
6.09e-03
7.66e-02
2.97e-01
6.69e-05
Obstacles
RNE
0.5061
0.6120
0.5410
0.6069
1.1956
0.0211
MSE
2.48e+03
4.98c+03
3.88c+03
1.95e+03
6.78c+03
1.76e+01
Laser Heat
RNE
0.1208
0.1711
0.1510
0.1071
0.1998
0.0102
4.2.1 CLASSICAL PDES
To validate the accuracy of DGNet as a general-purpose PDE solver, we first evaluate it on three well-established classical systems: the diffusiio--ominntted Allen-Cahn equation, the advection- driven Fisher-KPP equation, and the coupled FizzHug--Nauum system representing excitable me- dia. These tasks span distinct physical regimes (diffusion, advection-reaction. and multi-variable coupling), and are widely used as canonical benchmarks in scientific machine learning.
Table reports the quantitative results. DGNet consistently outperforms all baselines by large margin across all equations, achieving up to one to two orders of magnitude lower mean squared error (MSE). This indicates that the proposed discrete Green formulation, combined with physics- informed operators, substantially reduces sample complexity and delivers markedly improved accu- racy under limited supervision. Figure 3 provides qualitative comparisons. For the Allen-Cahi and Fisher-KPP equations, DGNet's predictions closely match the ground truth, whereas other models exhibit noticeable deviations, with MGN performing relatively better than the rest. The FitzHugh- Nagumo system is particularly challenging due to the formation of nonlinear spiral waves. Here, only DGNet successfully reproduces the propagation spiral structures over long horizons, while baselines either dissipate the patterns or introduce severe distortions. These results show that in a data-constr setting, DGNet can accurately solve PDEs with fundamentally different under- lying dynamics, highlighting its versatility and reliability data-efficient spatiotemporal solver.
8

---

Published as a conference paper at ICLR 2026
Ground truth
DGNet
DeepONet
MGN
MIPDE
BENO
PhyMPGN
Figure 3: Visualization of prediction results on classical PDE scenarios. Rows from top to bottom correspond to the Allen-Cahn, FitzHugh-Nagumo, and Fisher-KPP equations, respectively.
4.2.2 COMPLEX GEOMETRIC DOMAINS
We next examine the ability of DGNet to handle PDEs defined on irregular meshes with com- plex boundaries. These scenarios simulate contaminant transport in channel flow, with obstacles of varying shapes that induce rich flow structures. The three cases considered are: (i) a circu- lar cylinder, (ii) sediment deposits the walls, and (iii) combination of elliptical obstacles and airfoil-shaped wall structures. Such settings are challenging due to intricate boundary conditions and highly nonuniform meshes.
Table shows that DGNet achieves the lowest errors across all three domains, substantially out- performing baseline models. Figure provides qualitative comparisons. In the cylinder scenario, DGNet accurately captures the Kármán vortex street in the wake of the obstacle, while other mod- els exhibit distorted or dissipated vortical patterns. In the more complex sediment and obstacle cases, DGNet's predictions remain visually indistinguishable from ground truth, successfully re- producing contaminant filaments stretched and folded by the flow. By contrast, baseline methods either blur these fine-scale structures fail to track their spatial location. These results demonstrate that DGNet, benefiting from the strong structural prior provided by the discrete Green formulation, achieves data-efficient learning and robust performance even irregular geometries and nontrivial boundary conditions.
Ground truth
DGNet
DeepONet
MGN
MPPDE
BENO
PhyMPGN
Visualization of prediction results on scenarios with complex geometric domains. Rows from top to bottom correspond to the cylinder, sediment, and complex obstacle cases, respectively.
Figure 4:
4.2.3 GENERALIZATION TO UNSEEN SOURCE TERMS
salient property of DGNet is its ability to generalize beyond source terms encountered during training. To test this, we consider laser heat treatment task governed by the transient heat equation, where the source term f(x, t) corresponds to the spatiotemporally varying power density moving laser beams. During training, the model is exposed to trajectories generated from subset of source patterns. At test time, must adapt to entirely novel laser paths, such as spline and Lissajous curves, that were never capability to new forcing conditions, and serves as a stringent stress test data efficiency under limited training coverage.
during training. This setup directly evaluates the zero-shot generalization

---

Published as a conference paper at ICLR 2026
Quantitative results in Table show that most baseline models suffer severe degradation when fronted with unseen source terms: their errors increase by several orders of magnitude compared to the training distribution. In contrast, DGNet maintains stable accuracy, with almost performance drop relative to its results on seen source terms. Figure 5 further highlights this difference. DGNet accurately reproduces the spatiotemporal evolution of the temperature field, including the movement and spread of high-temperature regions, even under entirely novel laser trajectories. By comparison, baseline predictions collapse: the high-temperature regions become misplaced overly diffused, failing to follow the true laser paths. These experiments confirm that DGNet uniquely addresses critical limitation of existing neural PDE solvers, achieving robust generalization to unseen source terms- scenario ubiquitous in real-world scientific and engineering applications.
Ground truth
DGNet
DeepONet
MGN
MPPDE
BENO
PhyMPGN
Figure 5: Visualization of generalization performance on the laser heat scenario with unseen sources.
4.3 ABLATION STUDY
We compare DGNet against four variants on the complex obstacle scenario: (A) w/o Lphysisi re- moves the physics prior operator; (B) w/o Lneurat: removes the neural correction operator; (C) w/o Residual GNN: removes the residual GNN while retaining the discrete Green solver; (D) w/o Green: replaces the discrete Green solver with a generic end-to-end GNN that only consumes oper- ator features. The results are shown in Figure 6 Among all variants, w/o Green suffers the largest performance drop, highlighting the central role of the discrete Green solver as strong structural prior and inductive bias in the model. Removing Lphysies degrades performance more severely than removing Lneural. demonstrating that the physics prior provides essential structural knowledge while the neural correction mainly fine-tunes discretization errors. Finally, the decline in the w/o Resid- ual GNN variant confirms the residual GNN's effectiveness in capturing additional dynamics and improving overall accuracy.
0,8067
0.8
10-2
8.12--00
MISE
RNE
0.6-
100004
0.4
0.3198
1.49614
0.2
10
8.09 05
DGNN
w/f Lyphyrs
Resshui wijo Choen
DGNNtt
Lurnn w/lo Residual w/o Croen
Figure 6: Ablation study the complex obstacle scenario.
CONCLUSION
We introduced DGNet, a discrete Green network designed improve the data efficiency of neural PDE solvers through structural inductive bias. By discretizing Green's function on graphs, DGNet transforms a classical continuous tool into a practical learning framework that preserves the super- position principle and reduces the burden of rediscovering solution operator structure from data. Coupled with physics-based operators and GNN corrections, this architecture provides a princi- pled interface between physical fidelity and learnable flexibility. Extensive experiments across three categories of spatiotemporal PDE systems demonstrate that DGNet consistently achieves state-of- the-art accuracy under limited training trajectories, and maintains stable performance when general- izing to unseen source terms. Looking ahead, several open challenges remain, including extending the framework quasilinear PDEs where the superposition principle no longer holds, and scaling to large dimensionl systems. We view these directions as promising avenues for advancing discrete Green methods and broadening the impact in scientific machine learning.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENT
Q Yao supported by Beijing Natral Science Foundation (under Grant No. 4242039). Y. Wang sponsored by Beijing Nova Program.
REFERENCES
Keiiti Aki and Paul G Richards. Quantitative seismology. 2002.
John D Anderson. Computational fluid dynamics: the basics with applications. McGraw-Hill New York, 2002.
George Arfken, Hans Weber, and Frank E. Harris. Mathematical Methods for Physicists. Academic Press, 7th edition, 2013.
Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven dis- cretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116631):15344-1553 2019.
Suresh Bishnoi, Ravinder Bhattoo. Sayan Ranu, and N. M. Anoop Krishnan. Enhancing the induc- tive biases graph neural ODE for modeling physical systems. International Conference on Learning Representations, 2023.
Suresh Bishnoi, Jayadeva, Sayan Ranu, and N M Anoop Krishnan. BroGNet: Momentum- conserving graph neural stochastic differential equation for learning brownian dynamics. In In- ternational Conference on Learning Representations, 2024.
William E. Boyce. Richard DiPrima, and Douglas B. Meade. Elementary Differential Equations and Boundary Value Problems. Wiley, I 1tth edition, 2017.
Johannes Brandstetter, Daniel Worrall, and Max Welling, Message passing neural PDE solvers. In International Conference on Learning Representations, 2022.
Francisco Sahli Costabal, Simone Pezzuto, and Paris Perdikaris. A-PINNs: Physies-informed neu- ral networks on complex geometries. Engineering Applications of Artificial Intelligence, 127: 107324, 2024.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. In ICLR 2020 Workshop Integration of Deep Neural Models and Differential Equations, 2020.
Gwynne Evans, Jonathan Blackledge, and Peter Yardley. Numerical methods for partial differential equations. Springer Science & Business Media, 2012.
Mingquan Feng, Yixin Huang, Weixin Liao, Yuhong Liu, Yizhou Liu, and Junchi Yan. SINGER: Stochastic network graph evolving operator for high dimensional PDEs. In International Confer- ence on Learning Representations, 2025.
Robert Joseph George, Jiawei Zhao, Jean Kossaifi, Zongyi Li, and Anima Anandkumar: Incremental spatial and spectral learning of neural operators for solving large-scale PDEs. Transactions on Machine Learning Research, 2024.
John Guibas, Morteza Mardani, Zongyi Li, Andrew Tao, Anima Anandkumar, and Bryan Catanzaro. Efficient token mixing for transformers via adaptive fourier neural operators. In International Conference on Learning Representutions, 2021.
David Hahn and Necati Özisik. Heat conduction. John Wiley & Sons, 2012.
Charles Hirsch. Numerical computation of internal and external flows: The fundamentals of com- putational fluid dynamics. Elsevier, 2007.
Masanobu Horie and Naoto Mitsume. Physics-embedded neural networks: Graph neural PDE solvers with mixed boundary conditions. In Advances in Neural Information Processing Systems, pp. 23218-23222 2022.
11

---

Published as a conference paper at ICLR 2026
Tony Lelievre and Gabriel Stoltz. Partial differential equations and stochastic methods in molecular dynamics. Acta Numerica, 25:681-880, 2016.
Zhihao ei
et al. Harnessing scale and physics: multi-graph neural operator framework for PDEs on arbitrary geometries. In ACM SIGKDD Conference Knowledge Discovery and Data Mining, 2025.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar, Fourier neural oper- ator with learned deformations for PDEs general geometries. Journal of Machine Learning Research, 24(388):1 2023.
Rainald Löhner. Applied computational fluid dynamics techniques: An introduction based on finite element methods. John Wiley & Sons, 2008
Cheng Long, Siyuan Liu, and Zhiping Wang. PGODE: Towards high-quality system dynamics modeling. In AAAI Conference on Artificial Intelligence, volume 37, pp.
2023. Lu Lu, Pengzhan Jin, and George Em Karniadakis. Learning nonlinear operators via DeepONet.
Nature Machine Intelligence, 3(9):764-774, 2021.
Peter Lynch. The origins of computer weather prediction and climate modeling. Journal of Compu- tational Physics, 227(7):333133444 2008.
Xuhui Meng and George E Karniadakis. A composite neural network that learns from multi- fideli- data: Application to function approximation and inverse PDE problems. Journal of Computa- tional Physics, 401:109028, 2020.
Mark Meyer, Mathicu Desbrun, Peter Schröder, and Alan H Barr. Discrete differential -iirrnttillgcommetyy operators for triangulated 2-manifolds. In Visualization and Mathematics III, pp- 35-57. Springer, 2003.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzake and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021.
Stephen B Pope. Turbulent flows. Measurement Science and Technology, 12(11):2020-200 2021, 2001.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Roger Z Ríos-Mercado and Conrado portation systems: A statc-of-the-an review. Applied Energy, 147:536-555, 2015.
Optimization problems in natural gas trans-
Alvaro Sanchez-Gonzaler Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ODE integrators. In ICLR 2020 Workshop Integration of Deep Neural Models and Differential Equations, 2020a.
Alvaro Sanchez-Gonzalez Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning simulate complex physics with graph networks. In International Confer- 2020b.
ence on Machine Learning, pp.
Sungyong Seo, Chuizheng Meng, and Yan Liu. Physics-award difference graph networks for sparsely-observed dynamics. In International Conference Learning Representations, 2019.
Natarajan Sukumar and Ankit Srivastava. Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks. Computer Methods in Applied Mechanics and Engineering. 389:1 14333, 2022.
Allen Taflove, Susan C Hagness, and Melinda Piket-May. Computational electromagnetics: The finite-difference time-domain method. The Electrical Engineering Handbook, 3(629-670):15, 2005.
12

---

Published as a conference paper at ICLR 2026
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. In International Conference on Learning Representations, 2023.
Lloyd N Trefethen and Mark Embree. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, 2005.
Haixin Wang, Jiaxin Li, Anubhav Dwivedi, Kentaro Hara, and Tailin Wu. BENO: Boundary- embedded neural operators for elliptic PDEs. In International Conference on Learning Repre- sentations, 2024.
Gerhard Wanner and Ernst Hairer. Solving ordinary differential equations 11, volume 375. Springer Berlin Heidelberg New York, 1996.
Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: fast transformer solver for PDEs on general geometries. In International Conference on Machine Learning. pp. 53681-53705, 2024a.
Hao Wu, Huiyuan Wang, Kun Wang. Weiyan Wang. Yangyu Tao, Chong Chen, Xian-Sheng Hua, and Xiao Luo. Prometheus: Out-of-distribution fluid dynamics modeling with disentangled graph ODE. International Conference on Machine Learning, pp. 53870-53891, 2024b.
Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced physics-informed
neural networks for forward and inverse PDE problems. Computer Methods in Applied Mechanics and Engineering. 393:1 14823, 2022.
Bocheng Zeng, Qi Wang, Mengtao Yan, Yang Liu, Ruizhi Chengze, Yi Zhang, Hongsheng Liu, Zidong Wang, and Hao Sun. PhyMPGN: Physics-encoded message passing graph network for spatiotemporal PDE systems, In International Conference on Learning Representations, 2025.
A DERIVATION OF THE DISCRETE GREEN SOLVER
For completeness, we outline the derivation of the discrete Green solver and its relation implicit time integration. Unlike explicit propagation matrices, which are dense and costly to construct, our approach realizes the Green operator implicitly through sparse linear solves.
A.1 EXPLICIT vs. IMPLICIT GREEN FORMULATION
The discrete propagation matrix G(At) formally maps the state u to h+1 and serves as the discrete analogue of continuous Green's 's function:
u h+1 G(At)u..
Each column of G(At) corresponds the system response from a unit impulse at one node, while each row represents how contributions from all nodes aggregate to update a given node. However, explicitly forming GAtt) requires (N storage and computation, which infeasible for large meshes. We therefore derive an implicit realization based on stable time integration.
A.2 CRANK-NICOLSON DERIVATION
Starting from the semi-discretized ODE system
du Lu f,
dt
we apply the approximated by a central difference,
scheme for stability and second-order accuracy. The time derivative
du
u
t
dt
At
while Lu and f are approximated by their averages th and +++1 This yields
k+1
(Lu + ++1) + {(p + p*+1).
At
2
13

---

Published as a conference paper at ICLR 2026
Rearranging terms gives the sparse linear system
(I -I- -
(I- 2 ) (1* f
u²
which is equivalent to applying the discrete Green's 's function
G(At) = [I 1
A.3 EFFICIENT IMPLEMENTATION
In practice, do not compute the inverse explicitly. Since L is sparse and depends only the static mesh geometry, the coefficient matrix (I - L) is also sparse and remains fixed during rollout. We therefore adopt "factorize once, solve many times" strategy: single sparse LU factorization is performed before rollout, and hewarddbackwaar substitution reuses the factors at each time step. This reduces the per-step cost nearly linear in the number of nonzeros, making the discrete Green solver scalable to large systems.
B DATASET DETAILS
We evaluate DGNet on three categories of PDE systems: (1) classical equations (Allen-Cahn, Fisher-KIP KPP, FitzHug--Nagumo, (2) contaminant transport in complex geometric domains, and (3) laser heat treatment with varying source terms. High-fidelity simulation data are generated with the FEniCSx finite element library unstructured meshes and small time steps. Detailed parameters are reported in T
Below we outline the dataset settings.
B.1 CLASSICAL PDES
Allen-Cahn. Phase separation with 0.04 on [0, 1J.. A/ 0.005, 1000 steps. Mesh: perturbed 35 x 35 grid (~1296 nodes). Initial condition: random noise in [-0.5, 0.5]. Boundary condition: homogencous Neumann. Trajectories: 20 (10 train / 10 test).
KiPh Advection-reactie equation with p = 1.0, c [0.1, 0.12], domain J0, 1]3. At 0.00125, 1600 steps. Mesh: perturbed 40 x 40 grid (~1681 nodes). Initial condition: Gaussian blobs. Boundary: inflow Dirichlet, outflow natural. Trajectories: 20 (10 train 10 test).
FitzHugh-Nagumo. Two-component excitable system with parameters Da = 1/L.. D. 6.181/L3. r 9.657, = 0.5, B = 2.1, L, = 16x. Domain 10, At = 0.005, 1200 steps. Mesh: perturbed 99 x 99 grid (~10,0000 nodes). Initial: circular perturbation. Boundary: inflow Dirichlet (u = 0, 0), others Neumann. Trajectories: 10 (7 train / 3 test).
B.2 COMPLEX GEOMETRIC DOMAINS
All scenarios solve coupled Navier-Stoker and paaction-adveetion--dffiisssssss equations. Fluid den- sity P 1.0, viscosity Simulation: = 100s, At = 0.05s, 2000 steps. Trajectories: 20 (15 train / test).
5 10-3, contaminant /) 2 x 10 10-3 Pg 0.1, ka 0.025.
Cylinder. Rectangular channel [0, 24] [0, 8], obstacle: central circle (r 1.0), Re 240.
Sediments. Same channel, three asymmetric smooth bumps attached to top/bottom walls.
Complex Obstacles.
Same channel, two ellipses NACA airfoil-shapp shaped bumps/grooves.
B.3 LASER HEAT TREATMENT ( (VARYING SOURCE TERMS)
Transient heat conduction with = 7850, cp = 450, & = 50, Rgoev = 25. Domain: complex plate with gear cireular holes. Boundary: convective, initial T 298.15K. Heat source Q(z, t): 10 moving lasers with random paths (orbit, spline, Lissajous, etc.). Simulation: = 60s, At = 0.5s, 120 steps. Trajectories: 40 (20 train / 20 test).
14

---

Published as a conference paper at ICLR 2026
Table Simulation parameters of all datasets used in experiments.
Dataset
Domain
Nodes (approx.)
At
Steps
Traj. (train/test)
Classical PDEs
Allen-Cahn
[0, 1]2
1,296
0.005
1,000
10/
KiP
[0, 1
1,681
0.00125
1,600
10 10
Naguuo
[0, 1]
10,000
0.005
1,200
7/3
Complex Geometric Domains
Cylinder
0, 24]
10,8 8
2,275
0.05
2,000
15/5
Sediments
0, 24
x
5,758
0.05
2,000
15/5
Complex Obstacles
0, 24
x [0.8]
8,841
0.05
2,000
15/5
Varying Source Term
Laser Heat
Complex Plate
6,072
0.5
120
20 /0
C EXPERIMENTAL SETUP
All experiments conducted workstation with four NVIDIA RTX 4090 (24GB) GPUs. Train- ing typically completes within hours to one day depending on dataset size. Ground-truth PDE data are generated using FEniCSx (dolfinx 0.9.0). DGNet implemented in PyTorch 2.5.1 with CUDA 12.4 and Python 3.13.5.
C.I TRAINING HYPERPARAMETERS
All models are trained with Adam and step LR scheduler. We adopt trajectory segmentation for efficiency, and apply the loss only the first and last predictions of each segment to balance single-step accuracy and long-term stability. Table[4lists the hyperparameters for each experimental scenario (the three contaminant transport cases share identical settings).
Table 4: Training hyperparameters for each experimental scenario.
Scenario
Epochs
Decay Step
Decay
Sub-seq. Length
Train/Test Trajs.
Cahn
10/10
Fisher-KPP
24
Se-4
300
60
0.1
10
10/10
FitzHugh Nagumo
10
le-3
400
100
0.1
7/3
Cylinder
6
Se-4
300
100
0.2
15/5
Sediments
6
5o-4
300
100
0.2
15/5
Complex Obstacles
6
5e-4
300
100
0.2
15/5
Laser Heat
8
400
0.1
20/20
C.2 IMPLEMENTATION DETAILS
To solve the sparse linear system Eq. 6 Au inter
b, we employ a factorize-once strategy with
acceercelee sparse algebra.
Pre-computation. The system matrix A = (I- At L) depends only on mesh geometry and At, so it is built once before training. It is converted to CuPy sparse matrix and factorized using cupy. sparse. linalg. splu, producing LU factors cached as A.u..
Forward pass. At each step, the cached LU factors are used to efficiently solve A+11 by forward/backward substitution, avoiding explicit matrix inversion and enabling fast long-term rollouts.
Backward pass. We implement custom torch. autograd. Function with an adjoint for- mulation. Gradients are obtained by solving T go = gu,inter using the cached LU factors, ensuring differentiability and efficiency.
15

---

Published as a conference paper at ICLR 2026
D MORE EXPERIMENTAL RESULTS
D.1 VISUALIZATION OF LEARNED OPERATOR
A00
1.13
3000
0250
138
5.
-0005
Figure 7: Comparison of the learned operator with the reference. Left: input physical field. Center: reference operator computed by finite difference. Right: operator predicted by DGNet.
Figure compares the operator predicted by our model (right) with the reference computed by finite difference scheme (center), given the same input physical field (left). The learned operator closely matches the reference in both structural patterns and numerical scale, demonstrating that DGNet accurately captures operator-level dynamics.
Figure 8: Heatmap visualization of operator matrices. Left: corrected operator L.
Lphysics Middle:
Lncurat-
Right: final
Figure 8 shows heatmaps of Lphysics (left), Lncural (middle), and the final corrected operator L (right). The magnitude of entries in Lueural is significantly smaller than in Lphysies indicating that the neural component mainly provides fine-grained corrections top of the physical prior.
16

---

Published as a conference paper at ICLR 2026
D.2 INFERENCE COMPUTATIONAL OVERHEAD AND THROUGHPUT
We evaluate the computational efficiency of DGNet on single NVIDIA RTX 4090 GPU. Table 5 details the computational costs across all datasets. Notably, once factorized, the per-step inference highly efficient, achieving throughputs exceeding 90-130 steps/s across all scenarios.
Table 5: DGNet inference complexity and throughput measurements.
Dataset
Nodes
Pre-comp.
Inference
Total Inf.
Throughput Peak Mem.
(N)
Time (s)
Time (ms/step)
Time (s)
(steps/s)
(MB)
Allen-Cahn
1,296
0.6274
7.7419
8,3693
129.17
152.60
Fisher-KPP
1,681
3.6413
7.5764
15.7636
131.99
227.97
-aguum
10,000
10.7362
10.8723
23.7829
91.98
3830.72
Cylinder 2,275
2.5194
7.7895
18.0984
128.38
368.81
Sediments
5,758
6,0459
7.7943
21.6345
128.30
1570.49
Complex Obstacles
8,841
9.6110
9.7367
29.0844
102.70
3258.10 10
LaserHeat
6,072
1.9382
7.7363
2.8666
129.26
1570.65
To assess scalability, Figure illustrates the inference time and peak memory usage as the mesh size increases up to our largest experimental settings. The results demonstrate that DGNet maintains predictable and effective scaling behavior.
Scalability of DGNet with Increasing Mesh Size
20.0
5000
Inference Time
Peak Memory
4000
Inference Time (ms /stepp
15.0
Peak Memory (MB)
12.5
3300
10.0
2000
5.0
1000
0.0
2000
4000
6000
8000
10000
Number Nodes (N)
Figure 9: Scalability of inference time and peak memory with respect to mesh size.
Finally, Table 6 compares the single-trajectory inference time of DGNet against baseline models. Despite using pre-factorization, DGNet's inference time remains competitive with baseline meth- ods. These results confirm that DGNet achieves superior accuracy without imposing a significant computational burden compared existing methods.
17

---

Published as a conference paper at ICLR 2026
Table 6: Comparison of single-trajectory inference time (in seconds) between DGNet and baseline methods.
Dataset
DeepONett
MGN
MP-PDE
BENO
PhyMPGN
DGNet
Allen-Cahn
1.33
5.12
5.96
14.06
6.71
8.37
Fisher-KPP
2.77
12.74
19.20
75.92
17.40
15.76
FitzHuugh-aguume
106.52
15,48
45.23
67.32
10.72
23.78
Cylinder
6.38
20.15
24.90
133.01
12.47
18.10
Sediments
58.94
21.47
31.25
101.10
14.59
21.63
Complex Obstacles
141.92
23.59
47.36
108.90 9
20,31
29.08
Laser Heat
4.22
2.56
3.95
2.30
1.73
2.87
D.3 EXTENDED ABLATION STUDIES
To verify the robustness of our architectural components across diverse physical regimes, we extend the ablation study two additional scenarios: the FizZHug-NNauumo system and the Laser Heat experiment. We compare DGNet against the same four variants defined in Section 4.3 w/o Lneurai, w/o Residual GNN, and w/o Green GNN.
w/o Lphysics
The results, presented in Figure 10 and Figure are consistent with the findings in the main text. In both scenarios, the w/o Green GNN variant exhibits the most severe performance degradation, confirming that the discrete Green's function formulation is the fundamental backhone of our solver. Furthermore, removing the physics prior (w/o Ipphssica) leads a substantial increase in error, highlighting the necessity of embedding physical laws directly into the operator. Finally, both the neural correction (Lneurat) and the residual module contribute to the final accuracy, ensuring the model captures complex dynamics that pure discretization operators cannot resolve.
3.0
2. 2,704
2.5
7,24-06
MSD
10
6:204477
4.646-177
1.0
0U9471
1.18e
0.0662
10
DGNe
w/o Lypeee w/o Iuaan Wio Rooklul w/o Goen
DONet
w/o Iysas
w/lo braas
w/oo Risiinal w/o Gromm CNN
GNN
Figure 10: Ablation study results on the HarHHughNNagume scenario.
1.33400
0.20
0.1924
0,1539
0.15
M55
10P
5.90e++22
0.10
10F
4.00 01
0.05
w/o Lysyio
W/l Leaer w/o Resdual w/o Green
w/o Lptiis w/o L.vect M/oo Resdual w/o Green
CNN
GNN
Figure 11: Ablation study results on the Laser Heat scenario.
18

---

Published as a conference paper at ICLR 2026
D.4 EMPIRICAL SPECTRUM PLOTS
To empirically corroborate the stability analysis presented in Appendix we visualize the eigen- value spectra the operators on the FitzHugh-Naguma and Laser Heat datasets. In Figures 12 and 13 the blue markers represent the eigenvalues of the fixed physics prior Lpyyices while the orange markers represent the eigenvalues of the final learned hybrid operator
Lphgssics + Lacural.
The plots demonstrate that the cigenvalues of the total operator remain strictly confined to the left half of the complex plane (Re(A) < 0), respecting the stability boundary indicated by the red dashed line. This confirms that the neural correction neural functions as a stable perturbation, preserving the dissipative structure of the physical system while refining its dynamics, in alignment with our theoretical stability guarantees.
Eigenvalu Spectrum Operators (FHN)
threued
moogee Port oo Elppeaasss (1103)
Read Part Digenaaaae (Px3))
Figure 12: Eigenvalue spectrum analysis for the FitzHugh-Nagumo task. Blue indicates the physics prior physsccs and orange indicates the total hybrid operator L.
Eigenvalue Spectrum Operators (Heat Eq.)
6009
of Sprrstassss
0.02
0.00
0000
Bubere (Ragpe))
Par Elgeratin
Figure 13: Eigenvalue spectrum analysis for the Laser Heat task. Blue indicates the physics prior Lphysice and orange indicates the total hybrid operator L.
19

---

Published as a conference paper at ICLR 2026
D.5 SPEED-CCCURCY TRADE-OFF ANALYSIS
To further contextualize the efficiency of approach, analyze the speed-accuracy trade-off on the solver, which employs the Finite Element Method (FEM) with Newton nonlinear solver. By vary- ing the mesh resolution, we establish a baseline curve linking computational cost to simulation error.
dataset. We compare DGNet against a classical figh--iillity numerical
As shown in Figure 14 the classical numerical method exhibits steep increase in inference time as the error tolerance decreases. In contrast, DGNet demonstrates a characteristic "flat" cost curve; its inference time is determined by the network architecture and remains constant regardless of the achieved precision. Crucially, in the high-precision regime (MSE R 10-9), DGNet achieves comparable accuracy to the fine-mesh numerical solver but at a fraction of the computational cost. This highlights the distinct ad advantage of neural solvers: they do not aim to replace classical methods for unfnite-prreiitee tasks but serve as highly efficient surrogates within a practical accuracy range.
250
Numerical Method
DGNet (Ours)
200
Inference Time (s)
150
100
50
2
4
6
Mean Squared Error (MSE)
x10-6
Figure 14: comparison between DGNet and the classical FEM numerical solver on the FitzHugh-Nagumo task.
E THEORETICAL ANALYSIS OF NEURAL CORRECTION OPERATOR LNEURAL
This section focuses on the theoretical analysis of the learned neural correction operator Lncural- We highlight two core properties of this learnable component: Consistency and Stability. For the former, we show that the optimization process learns to cancel the local truncation error of the physical operator, thereby recovering the consistency of the discrete scheme. For the latter, establish spectral constraints for Lacural based on the Numerical Range. As long as the correction term satisfies specific norm bounds, the hybrid operator maintains the A-stability of the physical system.
We first provide clear definitions and the discrete scheme, similar Section Let Ω Rdd be bounded domain, and T7, be an irregular mesh with characteristic size h. Our graph is constructed based on this mesh Tr (using
subscript h to emphasize the characteristic scale). Discretizing the PDE 8,u- = Lu + ʃ into semi-discred ODE system on finite-dimensional space N the evolution of DGNet is governed by the hybrid operator Tap RNNNN
L Lphysics + Lncural (8).
where Lphysies is the fixed discretization operator based on mesh geometry. and (ne is the GNN correction term controlled by parameters 0. The single-step inference of the system follows
20

---

Published as a conference paper at ICLR 2026
the Srank-Ncoossn implicit update scheme:
At
(+ nº S())
k+1 1
In
T (++
2
where (AI) = (1- source term response.
At 2 lm)
-1 is the discrete Green's propagation operator, and S(f) is the
E.1 CONSISTENCY ANALYSIS
Our consistency analysis aligns with data-driven discretization framework proposed by Bar-Sinai et al. (2019). They demonstrated that learning correction coefficients on coarse grid effectively captures the unresolved sub-grid physics (truncation errors), enabling figh-fidelity simulation at low resolution.
Definition 1 (Spatial Truncation Error). Let Ps. be the projection operator from the continuous space to the discrete space. The local truncation error Tu( of the physical operator Lphysies with respect to the true operator C is defined as:
(u(u) Lphysics Pr(CC)).
On irregular meshes, due to mesh skewness, there is typically a non-zero low-order consistency error, i.e., ||Ta||| A O.
Theorem (Consistency Guarantee). For neural solver with universal approximation properties, if the training objective is to minimize the single-step prediction error J(0) k+1 + ||
(where denotes the model prediction and uGT denotes the ground truth), then the optimal neural
correction term 7 T satisfies:
neurall PAu = (u),
thereby eliminating the geometrically induced error of the physical discretization and recovering consistency.
Proof. Consider the single-step evolution from tu to 6x+1. The true solution uGT satisfies the semi- discrete equation with truncation error:
k+1
p ( + uGr) + 2 5 (u) +...
At
2
The predicted solution a
k++ of DGNet follows the update rule:
wer 1
2/2 (Uphysics + ( k+1 w + +...
At
Comparing the two equations, the consistency residual Econs of the hybrid operator relative to the true dynamics is:
Econs || Lpphysis + Lineraa) Pru- PnCull
[LLphysics PnCu) + Lncural
[TTr(u)= + (neuralP..ul]
During training, minimize the prediction loss [()), which is equivalent to finding correction operator such that ( ( + Ty (u). This implies that the neural network is explicitly driven learn the negative of the truncation error, thereby canceling the geometric error and recov- ering the consistency of the numerical scheme.
E.2 STABILITY ANALYSIS
Our stability analysis relies the theory of Numerical Range. Following Trefethen & Embree (200)), we use the spectral enclosure property of the numerical range to bound the system's dy- namics. To establish strict energy stability, we employ the Logarithmic Norm of the operator. As detailed in Wanner & Hairer (1996), negative logarithmic norm ju 4n) < 0 is a sufficient condition for unconditional stability in stiff systems. Our proof based on this fact.
Definition 2 (Numerical Range and Logarithmic Norm). The numerical range of a matrix A is (x.Az) C \ {0)}, which possesses the following key properties:
defined as W (A)
(0,0)
21

---

Published as a conference paper at ICLR 2026
Spectral Enclosure: All eigenvalues X(A) of the matrix are contained within the numerical range, Lee, (A) S (A).
Dissipativity Bound: The supremum of the real part of the numerical range equals the logarithmic norm: sup (o(W(A)) (2) Amax A+44
2
Theorem (Stability Guarantee). Assume the physical operator Lphysies is coercive (dissipative), i.e. for any x + O, Re(r, Lphysiicat) S -n||>r|| (where vj > determined by the physical system). If the output of the neural correction operator Lmnural is constrained by architecture (e.g.. bounded activation function tanh) such that its spectral norm satisfies L S 7, and the condition <
holds, then the hybrid operator Ln Lpbysics + Licurat satisfies:
Spectral Stability: The real parts of all eigenvalues are strictly negative, i.c., Re(A(Ln)) < O.
Energy Stability: The discrete evolution operator (I - - (1 At Ln) is strictly contrac- At tive, ie., || At Ln) -1 (I+ 2 Ln) 112 <.
Proof. For any unit vector x <1<0112 = 1), the quadratic form corresponding to the hybrid operator is:
(or, Lnx) (xr, Lphysiesat) + (ar, L)eara
Taking the real part and using the Cauchy-Schin Schwarz inequality:
Re(æ, Lnx) Re(x, Lphysiess) + Re(x, L)curat
Due to the coercive dissipativity of Lppysiic and the boundedness of Lneerall
Re(x, Lax) I -7 + (rx, [meeratt:| S -1 + | S -1) +
This indicates that the numerical range (L,) lies entirely to the left of the line x = -n + y in the complex plane. Therefore, by the spectral enclosure property (Ln)) S (Ln),, we directly obtain the upper bound for the real parts of all eigenvalues:
max Re(A) s sup Re(W - (Ea(Ia)
This proves the spectral stability of the operator TTp
For energy stability, note that the norm of the Crank-Nicolson propagation operator (ل (I 2 Ln) satisfies the following bound:
(T
Al In
-1
Al
1+ At u2 (Ln)
1
In
1l2
2
1 Al (2()
Since (2())) sup Re(W < -1) + < 0, we know that the above norm is less than 1, obtaining energy stability.
Therefore, as long as the magnitude of the neural correction term is smaller than the physical dis- sipation rate 77, the hybrid operator of DGNet possesses a strictly left-half sat-hall-plane spectral distribution and unconditional A-stability.
22

---

Published as a conference paper at ICLR 2026
F SUPPLEMENTARY EXPERIMENT
EXPERIMENTS WITH 3D SCENES
To evaluate DGNet's scalability to high-dimensionis geometries and complex boundary conditions, constructed 3D transient heat conduction benchmark.
Geometry and Physical Model. The domain is a mechanical bracket (L-shaped support with a central bore) discretized into unstructured tetrahedral mesh with approximately 6,000 nodes (see Figure 15). We solve the transient heat equation with convective cooling:
pep (kVT) hcony (T Tamb) + Q(xx(.)).
-3, Cp 450 J -l- KK_1 50 W m - and using the following parameters: 7850 kg m-
hoony
25 W m -K--1 The simulation spans ' 60,s with a time step At 0.5,s.
06
0.8
02
0.6
0.2
Figure 15: Visualization of the 3D heat conduction scenario, showing the mechanical bracket geom- etry.
Laser Source Model. The source Q t) models a superposition of 10 moving Gaussian laser spots. To rigorously test generalization, we generate trajectories that combine deterministic paths (mimicking industrial processing) and stochastic paths (random splines and reflections). The beam radius is randomized in [0.5, 2.5).mm. mm, and power profiles are temporally ramped to create multi-modal forcing.
We generated 60 trajectories in total, split into 20 for training and 40 for testing. The experimental results are shown in Table 7 DGNet significantly outperforms baseline methods. This confirms that framework generalizes effectively to complex 3D geometries.
Dataset and Results.
Table 7: Quantitative results the 3D laser-heat task. DGNet achieves the lowest MSE and RNE compared to baselines.
Metric
DeepONet
MGN
MP-PDE
BENO
PhyMPGN
DGNet
MSE
4.4le+03
4.60e+02 0.1249
2.59e+04 0.1993
5.28e+02 0.0938
8.83e+03 0.2352
4.42e+01 0.0314
RNE
0.1630
23

---

Published as a conference paper at ICLR 2026
F.2 EXPERIMENTS WITH COMPLEX GEOMETRY AND HIGHER REYNOLDS NUMBERS
To further test the robustness of the model, we constructed more challenging variant of the Con- taminant Transport scenarios presented in Section 4.2.2 Experimental Setup. We maintained the identical governing equations, source injection strategy, and time-step settings as the previous ex- periments. The modifications were strictly limited to aspects:
Complex Boundaries:
The channel walls were modified from straight to sinusoidal (wavy) forms force continuous flow redirection.
Higher Reynolds Number:
The Reynolds number increased to Re 300, pushing the flow further into the transitional regime with richer vortex shedding dynamics.
Figure 16: Visualization of the modified wavy channel scenario (Re 300).
Dataset and Results. We generated trajectories under this configuration, and the scene visual- ization is shown in Figure 16 Experimental results shown in Table 8 Despite the more complex geometry, DGNet maintained robust performance.
Table 8: Quantitative comparison on the High-Reynolds (Re 300) fluid and pollutant transport task.
Metric
DeepONet
MGN
MP-PDE
BENO
PhyMPGN
DGNet
MSE
2.38c-01 0.8539
3.59c-03
7.94c-03 0.4729
2.30c-02 0.6292
7.95c-02 1.1838
6.29e-04 0.0385
RNE
0.3591
F.3 EXPERIMENTS WITH LARGE-SCALE LIGHT-DRIVEN REACTIONS
To validate the model's s scalability and generalization capability, we introduce a Light-Driven Chem- ical Reaction benchmark.
Figure 17: Snapshot of the light-driven reaction scenario.
Simulation Setup. The domain is circular Petri dish (R = 1) discretized into a fine unstructured mesh containing approximately 30,000 nodes. We solve transient reaction-diffusion equation with zero-flux boundary conditions:
ac
DVPC &C + f(x, t),
0t
24

---

Published as a conference paper at ICLR 2026
where D = 5 10-4 and the decay rate kt = 5 10-3
Complex Dynamic Source. The source term f(x. t) represents dynamic light intensity driving the reaction. To prevent overfitting to simple distributions, six distinct dynamic patterns (including rotating spirals, pulsating rings, moving spots, and traveling waves) with randomized parameters. This results in a highly nonlinear and unpredictable spatiotem- poral field.
construct the source by superimposing
We significantly expanded the dataset size 200 trajectories (40 for training, 160 for testing), spanning 7 = 100s with Á/ = 0.5s. The experimental results are shown in Table DGNet demonstrates the best accuracy and scalability. It scale to large systems and generalize to unseen complex source terms.
Dataset and Results.
Table 9: Quantitative results on the large-scale light-driven reaction task.
Metric
DeepONet
MGN
MP-PDE
BENO
PhyMPGN
DGNet
MSE RNE
6.11e+03 0.8832
1.93++02 0.1582
1.37e+01 0.0422
6.04e+02 0.2796
2.87e+02 0.1921
3.97e+00
0.0227
F.4 MORE SOURCE TERM GENERALIZATION EXPERIMENTS
To evaluate DGNet's robustness to variations in source characteristics in zero-shot setting (the trained model is fixed), performed robustness tests that vary three key source attributes: spatial bandwidth (spot size), amplitude (peak power), and the number of concurrent sources. Each attribute was scaled relative the training distribution by factors {0.5, 1, 2, 3} and evaluation was performed on the Laser-Heat 3D scenario.
Table IO summarizes representative numeric results (MSE) under these perturbations. DGNet gen- eralizes effectively across the tested range. These findings indicate strong robustness to changes in source characteristics.
Table 10: Robustness to source characteristics. Each column scales the corresponding attribute by the indicated factor relative training.
Source Attribute
0.5
1
2x
3x
bandwidth
1.53e+01
1.76e+01
3.88e+01
7.36e+01
amplitude
2.27e+01
1.76e+01
1.62e+01
1,95e+01
number of concurrent sources
1.09c+01
1.76c+01
3.26c+01
4,40c+01
F.5
FAILURE MODE: QUASILINEAR PHASE TRANSITION
To probe regimes where the linear superposition assumption may break, we constructed a controlled failure-mode experiment based on a quasilinear phase-transition heat equation:
pcen (T),T V (k(T)VT) + Q(x,t).
(11)
where both the effective heat capacity Cefr(T) and thermal conductivity k(T') undergo sharp, temperature-depeedentt changes (e.g., step changes modeling phase transitions or abrupt material property shifts). This introduces strong nonlinear coupling between temperature and transport that violates linear superposition.
We employed temperature T.. Figure 18 presents visual comparison between the ground truth and the DGNet prediction. As observed in the Ground Truth (left), the temperature field exhibits a distinct plateau (uniform red region). This accurately reflects the physical mechanism phase transition, where the effective heat capacity (T')(TT) surges at T.,, causing the system to absorb input energy as latent heat rather than increasing the local temperature. In contrast, the DGNet Prediction (right) fails to capture this nonlinear saturation effect. Instead, it predicts diffusive distribution with high- temperature hotspots (bright yellow regions). This behavior occurs because the architecture relies on the superposition of source responses; it treats the energy injection additive linear heating.
ceteraatureddpendeeen parametrizations with a sharp transition centered near a critical
25

---

Published as a conference paper at ICLR 2026
Ground Truth
Prediction
Figure 18: Visualization of the failure mode in the quasilinear Phase Transition scenario.
ignoring the drastic, state-dependent changes in material properties. Consequently, prediction errors concentrate heavily around the phase transition interface and regions of source overlap.
These results empirically confirm boundary of our method: while DGNet is highly effective for semilinear systems, the linear superposition assumption limits its applicability in quasilinear regimes characterized by strong state-dependent parameter coupling.
F.6 COMPARISON WITH RECENT BASELINES
To comprehensively position DGNet against recent advancements, we evaluated three baselines SINGER (Feng et al. 2025), AMG-GNN (Li et al. 2025), and Transolver (Wu et al. (024a). which represent multi-scale architecture, graph- based architecture and Transformer architecture, respec- tively.
The results summarized in Table 11 show that DGNet achieves substantially lower error rates across all evaluated metrics, while recent baselines exhibit noticeably degraded performance under the limited-data setting considered in this work. As discussed earlier, we attribute this performance gap primarily to data efficiency. Although these data-driven models are powerful and expressive, they typically rely on thousands of training trajectories to accurately capture complex dynamics. In contrast, by embedding physical priors through the discrete Green's ''s function, DGNet is able to gen- eralize effectively from limited data, whereas purely data-driven baselines lack sufficient inductive bias.
Table 11: Quantitative comparison with recent state-of-the of-the-art baselines (SINGER, AMG-GNN, Transolver) representative tasks.
Dataset
Metrie
SINGER
AMG (Harnessing)
Transolver
DGNeet
MSE
4,58e-06
3,30e-05
7.87e-06
1.18c-07
FitzHughNaguumo
RNE
1.3844
2.4837
0.6729
0.0952
MSE
9.41e-03
1.40e-03
8.94e-04
6.69c-05
Complex Obstacles
RNE
0.7746
0.3058
0.1855
0.0211
MSE
9.58c+03
2.96e+03
5.72c403
1.76e+01
Laser Heat
RNE
0.2504
0.1831
0.1785
0.0102
26


Reviews:
[{"review_id": "L5h5a3beOm", "reviewer": "Reviewer_xgKd", "summary": "This paper proposes DGNet, a neural solver for spatiotemporal PDEs that explicitly decouples the system evolution from the source response by turning Green’s function theory into a discrete, graph-based update rule. Starting from a Crank–Nicolson semi-discretization, the authors derive a one-step update that can be interpreted as applying a discrete Green operator to propagate the current state and accumulate the effect of the source over the time slab. Experiments span classical PDEs, transport on irregular domains, and generalization to unseen source terms. DGNet reports SOTA accuracy across tasks and especially large gains when tested on entirely novel sources; an ablation suggests the discrete Green solver is the key driver of performance.", "questions": "1. Stability & constraints on $L_{\\text{neural}}$.\nCan you provide spectral/energy bounds or constraints on $L_{\\text{neural}}$ (e.g., skew-symmetry, diagonally dominant correction, Lipschitz bounds) that preserve the A-stability properties of the CN-like update when combined with $L_{\\text{physics}}$? A short lemma or empirical spectrum plots across datasets would strengthen confidence.\n\t\n2. Ablations beyond one geometry.\nYour ablation indicates “w/o Green” is worst on complex obstacles. Could you add two more scenarios (one classical PDE, one laser heat) to separate the contributions of $L_{\\text{physics}}$, $L_{\\text{neural}}$, and the residual GNN more generally? This would clarify the portability of each component. \n\t\n3. Runtime, memory, and scaling.\nPlease include a main-text table: per-step wall-clock, memory footprint, and speedup vs. baselines at equal accuracy, and a mesh-scaling plot up to your largest meshes. What happens if the sparse LU factorization is replaced with ILU/AMG preconditioning for larger systems?\n\t\n4. Robustness to source characteristics.\nThe laser task varies paths; could you report sensitivity to source bandwidth/amplitude and number of concurrent sources (e.g., 1 vs. 10 lasers), including failure modes where superposition might break due to strong nonlinearity?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Solid derivation of a discrete Green update tied to an implicit Crank–Nicolson step, with practical “factorize once, solve many” sparse solves; the operator depends only on mesh and \\Delta t, enabling reuse during rollout.\n\n- Well-motivated hybrid operator that anchors learning with physics priors (Green–Gauss gradient, cotangent Laplacian) plus GNN corrections; clear architectural overview.\n\n- Broad empirical coverage with consistent SOTA on both log-MSE and relative $\\ell_2$ across diverse regimes; the table highlights DGNet edges on every scenario.", "weaknesses": "- Some implementation details (e.g., GPU sparse LU and custom adjoint) are in the appendix; a brief complexity/throughput table in the main text would help readers assess practicality across datasets.\n\n- Uses standard message-passing blocks; novelty is concentrated in the Green discretization + operator split, not in GNN mechanics.  \n\n- Although there is an ablation, the residual GNN’s incremental benefit and the sensitivity to factorization accuracy (e.g., fill-in thresholds, preconditioners) could be quantified more rigorously beyond one scenario.\n\n- While the CN-derived update is principled, the scope of PDE classes is primarily parabolic/weakly nonlinear; truly strongly nonlinear or hyperbolic systems are deferred to future work.\n\n- Stability/consistency guarantees for the learned correction are not theoretically analyzed.\n\n- Benchmarks are 2D; without a 3D case, claims about broad scientific impact are somewhat aspirational.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "YYtHA4EKYG", "reviewer": "Reviewer_HxRT", "summary": "This paper proposes DGNet, a discrete Green network framework for solving spatiotemporal PDEs with a particular focus on explicit modeling of source terms.\nUnlike conventional neural PDE solvers that implicitly mix the source and system state, DGNet explicitly decouples the system evolution from the source response, inspired by the Green’s function formalism.\nThe model constructs a graph-based discrete operator that preserves the superposition principle. The operator $\\mathcal{L}$ is composed of a physics-based component $\\mathcal{L}{\\text{physics}}$ and a neural correction component $\\mathcal{L}{\\text{NN}}$, where the former computes spatial derivatives such as gradients and Laplacians numerically, and the latter leverages a message passing neural network (MPNN) to correct approximation errors.\nDGNet further incorporates Crank–Nicolson time integration to ensure stability in temporal evolution. Across multiple PDE benchmarks, including irregular mesh scenarios and novel source terms, DGNet achieves state-of-the-art accuracy and demonstrates strong robustness compared to existing operator-learning methods.", "questions": "- In the ablation study, does “w/o $\\mathcal{L}_{\\text{NN}}$” correspond to solving the PDE purely via numerical methods (i.e., without any learned correction)?\n- How exactly is $\\mathcal{L}_{\\text{physics}}$ computed in practice—are the gradients and Laplacians incorporated as feature inputs to the MPNN, or directly as numerical differential operators?\n- The paper emphasizes strong generalization to unseen source terms, but given the small test sizes (e.g., 3–20), can the authors provide additional experiments or statistical evidence to support this claim?\n- How does DGNet compare to PHYMPGN (ICLR 2025) under identical experimental settings? The differences in reported performance seem unexpectedly large.\n- Could larger-scale experiments (more trajectories or longer temporal horizons) be conducted to further validate generalization?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper addresses a longstanding limitation in neural PDE solvers—their inability to generalize to unseen source terms—by explicitly incorporating the Green’s function concept into a learnable, graph-based framework.\nThis formulation is conceptually elegant: by separating the effect of the source term from the system dynamics, DGNet captures the response structure in a principled and interpretable manner.\nThe combination of physics-based discretization and neural correction strikes a strong balance between numerical fidelity and data-driven flexibility. In particular, the hybrid operator $\\mathcal{L} = \\mathcal{L}{\\text{physics}} + \\mathcal{L}{\\text{NN}}$ effectively merges computational stability with adaptability, while maintaining computational efficiency.\nEmpirical results show exceptionally high performance—often outperforming existing baselines by large margins—and stability on unseen forcing terms.\nOverall, DGNet’s design provides a conceptually grounded and practically effective approach for learning PDE dynamics on irregular meshes and under novel source conditions.", "weaknesses": "Despite its strong results, the paper suffers from several clarity and interpretability issues that limit its accessibility.\nFirst, the presentation of the operator $\\mathcal{L}$ and its components lacks precision. Although equations define $\\mathcal{L}{\\text{physics}}$ and $\\mathcal{L}{\\text{NN}}$, it remains unclear how $\\mathcal{L}{\\text{physics}}$ is numerically computed and integrated into the overall update rule—whether gradients and Laplacians are used merely as input features or as direct numerical operators.\nSimilarly, the interaction between $\\mathcal{L}{\\text{NN}}$ (the MPNN correction) and the source term $f$ is only superficially discussed, leaving readers uncertain about how the model actually combines physical and learned dynamics to advance the PDE state $u$.\nAlthough the abstract emphasizes source-term generalization, the main text provides limited explanation of how the method ensures this capability beyond architectural intuition.\nSecond, the neural component is relatively minimal. Apart from the residual MPNN correction, the rest of the solver heavily relies on standard numerical discretizations. While this hybrid structure is defensible, the role of the NN component could be elaborated to justify the “learning” aspect of the framework.\nThird, the evaluation protocol raises concerns about robustness. Reported test sets are extremely small (3–20 samples), making it difficult to confidently assess generalization claims. The performance gains may partially stem from limited sampling rather than true model generalization.\nFinally, after the numerical solution step, the exposition becomes sparse, with several mathematical derivations (Eqs. 5–7) presented without intuitive interpretation or discussion.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "R3zXDkzL5a", "reviewer": "Reviewer_dKtN", "summary": "The work proposes DGNet, a GNN-based model, to solve a class of PDEs. The method provides an explicit way to decouple the system evolution term from the effect of the source term based on the superposition principle of the PDE theory. The authors, in addition, investigate the relation between the Green’s function and a discrete time integration method, yielding the discrete counterpart of the superposition using the Green’s function. Further, the spatial derivative operator is expressed as a mixture of the classical numerical method and GNN. The experimental results show that the proposed method demonstrates the best accuracy among the considered baselines and generalization to unseen source terms. The authors conduct an ablation study to confirm that the proposed ingredients work properly.", "questions": "1. The method assumes $u$ represents a scalar field (line 125), but can we extend it to vector and tensor fields?\n2. How is the source term varied in the datasets? Can the method generalize to a source with a norm larger than that of the training dataset (e.g., a two-times larger source)? If not, the authors should explicitly present the limitations and the cases where the proposed model fails.\n3. The work assumes the mesh is generated using the Delaunay triangulation. Why did the authors choose that configuration? Can we consider using other types of meshes, e.g., quadrilateral and polygonal meshes? How about 3D meshes? What would be the difficulty in extending the method?\n4. The paper argues for the necessity of the neural correction operator for coarse or skewed meshes (line 265), but since the mesh is generated using the Delaunay triangulation, the reviewer considers there to be no skewed mesh in the present work. In addition, in the experiments, in particular, the FitzHugh–Nagumo dataset, the mesh resolution seems high enough. In that case, how does the author explain the necessity of the neural correction?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The focus on the source term is interesting. Since the source/forcing term is used across a wide range of research and engineering applications, the work is of great importance to the community.\n2. The method is simple and reasonable. The authors successfully leverage classical theory to develop an effective approach for elegantly handling the source term.", "weaknesses": "1. The applicable domain and limitations of the work are unclear. Since the method relies on the superposition and linear approximation of the spatial operator $\\mathcal{L}_\\boldsymbol{x}$, the reviewer considers it to work only for semilinear PDEs, although the paper says “PDE” with no specifications. The authors should clarify the assumptions and limitations of the work.\n2. Section 3.4.1 is unclearly written. The paper uses $A_i$ and $d_i$ for volume and area, respectively, but this is outside the standard notation. In addition, there is no explanation about $N(i)$, $\\alpha_{ij}$, and $\\beta_{ij}$. The reviewer recommends adding a figure showing the variable settings for improved readability.\n3. The experimental evaluation is weak. The work uses machine learning methods as baselines, but the comparison of computational speed and accuracy with a classical numerical solver should also be provided. In particular, since the method includes the scheme from the classical solver, it is not obvious that the proposed approach is more efficient than its classical counterpart. In addition, since machine learning methods have some error, the evaluation of a classical solver should be performed across multiple settings, e.g., varying the mesh resolution and convergence threshold, to capture the speed-accuracy tradeoff.\n\nMinor point:\n\n* Labels in Figure 5 and description in the main text are not aligned (e.g., w/o fit in the table vs w/o Residual GNN in the text). The reviewer recommends aligning them in the same word.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DS8AxKnwMq", "reviewer": "Reviewer_pbbn", "summary": "The paper proposes **DGNet (Discrete Green Network)** for learning **spatiotemporal PDEs** on irregular meshes, with a core design that **explicitly decouples system evolution from source response** via a **discrete Green’s function** formulation. Concretely, it discretizes ( $\\partial_t u = L_x[u] + f$ ) into a **Crank–Nicolson** update and realizes a per-step Green operator ($G(\\Delta t) = (I-\\tfrac{\\Delta t}{2}L)^{-1}$), enabling superposition of (i) state propagation and (ii) source-term response. The spatial operator is built as a **physics–neural hybrid** ($L=L_{\\text{physics}}+L_{\\text{neural}}$): ($L_{\\text{physics}}$) uses mesh-aware gradient/Laplacian discretizations, while ($L_{\\text{neural}}$) is a GNN correction for discretization errors; a lightweight residual GNN further captures leftover dynamics. An efficient **factorize-once, solve-many** sparse LU scheme makes long rollouts practical. Experiments across **classical PDEs**, **complex geometries**, and **unseen source terms** show consistent SOTA accuracy, with marked robustness when test-time sources differ from training.", "questions": "1. **Precomputation & fairness.**\n   Please report a cost table (params, FLOPs/step, GPU memory, LU pre-factorization time for Eq. (9), wall-clock/epoch, hardware) and confirm matched training budgets across baselines.\n\n2. **Scalability & factorization reuse.**\n   Results on substantially larger meshes (and ideally 3D) would align with the motivation. Also clarify when (L) changes across samples (coefficients/BCs/mesh): must you re-factor each time, and what is the end-to-end impact?\n\n3. **Benchmarks & baseline selection.**\n   Add stronger, recent multiscale/transformer/graph operator baselines for irregular spatiotemporal PDEs, and justify inclusion of elliptic-focused baselines. For **cylinder flow**, reconcile why baseline errors are far worse than commonly reported (detail resolution/Re/splits/targets/configs).\n\n4. **Robustness & hybrid ablations.**\n   Provide stress tests (larger $(\\Delta t)$, stiff sources, mesh noise/density shifts) and ablations separating $(L_{\\text{physics}})$ vs. $(L_{\\text{neural}})$ to show stability and the benefit of the physics–neural hybrid.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. **Addresses challenging spatiotemporal PDEs on irregular meshes.**\n   The setting combines *both* complex geometry (unstructured meshes) and temporal evolution, which is closer to real CFD/physical simulations than static or grid-regular cases. Tackling this regime is practically meaningful. \n\n2. **Green’s-function–motivated decomposition.**\n   By formulating updates through a discrete Green operator, the method **explicitly separates** (i) the evolution of the initial state and (ii) the accumulated response to time-varying sources. This superposition-friendly design improves interpretability, allows cleaner handling of unseen source patterns, and provides a principled bridge between physics and learning.\n\n3. **Physics-informed design on irregular meshes.**\n   The hybrid operator and loss terms embed **physically consistent priors** (e.g., mesh-aware differential operators, stability-friendly time-stepping), encouraging fidelity to the governing equations while letting the learned components correct discretization errors. This tends to enhance robustness, boundary handling, and generalization across meshing/sampling changes.", "weaknesses": "1. **Heavy reliance on pre-factorization; fairness not quantified.**\n   The method hinges on a **“factorize once, solve many”** sparse LU of ($(I-\\tfrac{\\Delta t}{2}L)$) before rollout, then reuses the factors each step. While efficient per step, this adds a non-negligible **precomputation and memory** burden uncommon in NN-only baselines, yet there’s **no cost table** (params/FLOPs/GPU memory/wall-clock) to establish fairness. The paper mentions hardware (4×RTX4090) and vague training durations (“hours to one day”) but lacks matched-budget reporting across methods. A standardized efficiency comparison is needed.  \n\n2. **Small problem sizes vs. stated motivation.**\n   Core 2D meshes are relatively modest (e.g., **Cylinder ~2.3k nodes**, Sediments ~5.8k, Complex Obstacles ~8.8k; classical PDEs 1.3k–10k). These scales undercut the claim of handling challenging irregular domains; results on **larger meshes/3D** or higher Reynolds/complex boundaries would better substantiate practical robustness. \n\n3. **Baseline scope/positioning leave doubts.**\n   The baseline set (DeepONet, MGN, MP-PDE, PhyMPGN, **BENO**) mixes operator, graph, and hybrid methods, but (i) it’s unclear **how comparability is enforced** beyond a sentence (“tuned to have comparable parameters/training costs”) and (ii) some choices may be **mismatched to spatiotemporal settings** (e.g., BENO targets elliptic PDEs), while **strong recent multiscale/graph/transformer/point-operator** [1-5] baselines on irregular domains are missing. A head-to-head with closer, stronger contemporaries and a **capacity-controlled** comparison would sharpen novelty and fairness. \n\n[1] Feng, Mingquan, et al. \"SINGER: Stochastic Network Graph Evolving Operator for High Dimensional PDEs.\" The Thirteenth International Conference on Learning Representations. 2025.\n\n[2] Zhang, Xuan, et al. \"SineNet: Learning Temporal Dynamics in Time-Dependent Partial Differential Equations.\" The Twelfth International Conference on Learning Representations.\n\n[3] Wang, Qi, et al. \"P $^ 2$ C $^ 2$ Net: PDE-Preserved Coarse Correction Network for efficient prediction of spatiotemporal dynamics.\" Advances in Neural Information Processing Systems 37 (2024): 68897-68925.\n\n[4] Li, Zhihao, et al. \"Harnessing scale and physics: A multi-graph neural operator framework for pdes on arbitrary geometries.\" Proceedings of the 31st ACM SIGKDD Conference on Knowledge Discovery and Data Mining V. 1. 2025.\n\n[5] Wu, Haixu, et al. \"Transolver: A Fast Transformer Solver for PDEs on General Geometries.\" International Conference on Machine Learning. PMLR, 2024.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)