cs.LGMar 25, 2026

Learning Mesh-Free Discrete Differential Operators with Self-Supervised Graph Neural Networks

Authors: Lucas Gerken Starepravo, Georgios Fourtakas, Steven Lind, Ajay B. Harish, Tianning Tang, Jack R. C. King

Organizations: School of Engineering, The University of Manchester, Manchester, UK · School of Engineering, Cardiff University, Cardiff, UK

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Mesh-free numerical methods provide flexible discretisations for complex geometries; however, classical meshless discrete differential operators typically trade low computational cost for limited accuracy or high accuracy for substantial per-stencil computation. We introduce a parametrised framework for learning mesh-free discrete differential operators using a graph neural network trained via polynomial moment constraints derived from truncated Taylor expansions. The model maps local stencils relative positions directly to discrete operator weights. The current work demonstrates that neural networks can learn classical polynomial consistency while retaining robustness to irregular neighbourhood geometry. The learned operators depend only on local geometry, are resolution-agnostic, and can be reused across particle configurations and governing equations. We evaluate the framework using standard numerical analysis diagnostics, showing improved accuracy over Smoothed Particle Hydrodynamics, and a favourable accuracy-cost trade-off relative to a representative high-order consistent mesh-free method in the moderate-accuracy regime. Applicability is demonstrated by solving the weakly compressible Navier-Stokes equations using the learned operators.

Figures & tables

Figure 1: We learn the mapping from relative position of particles within a local neighbourhood, to a local set of weights that approximate differential operators. The learned operator is physics-agnostic, and local to the computational stencil. The framework consists of a lifting the relative positions to a higher dimension with a neural network, followed by a stack of message-passing graph layers, and a final output head that maps latent representations to operator weights.
Figure 2: Second-order normalised learned operators, colours indicate weight value. Models were trained with particle disturbance of ϵ=1.0\epsilon=1.0 and inferred with the same noise.
OperatorMetricxxyyx2/2x^{2}/2xyxyy2/2y^{2}/2
NeMDOp=2x\text{NeMDO}_{p=2}^{x}MAE8.84×10−58.84\times 10^{-5}7.93×10−57.93\times 10^{-5}5.77×10−55.77\times 10^{-5}5.59×10−55.59\times 10^{-5}4.74×10−54.74\times 10^{-5}
St.d.1.14×10−41.14\times 10^{-4}1.04×10−41.04\times 10^{-4}7.38×10−57.38\times 10^{-5}7.35×10−57.35\times 10^{-5}6.25×10−56.25\times 10^{-5}
NeMDOp=2Δ\text{NeMDO}_{p=2}^{\Delta}MAE5.22×10−45.22\times 10^{-4}5.09×10−45.09\times 10^{-4}3.69×10−43.69\times 10^{-4}4.09×10−44.09\times 10^{-4}3.62×10−43.62\times 10^{-4}
St.d.6.91×10−46.91\times 10^{-4}6.69×10−46.69\times 10^{-4}4.94×10−44.94\times 10^{-4}5.49×10−45.49\times 10^{-4}4.77×10−44.77\times 10^{-4}
Wendland C2x\text{Wendland C2}^{x}MAE5.14×10−25.14\times 10^{-2}4.01×10−24.01\times 10^{-2}2.70×10−22.70\times 10^{-2}2.63×10−22.63\times 10^{-2}1.36×10−21.36\times 10^{-2}
St.d.6.21×10−26.21\times 10^{-2}5.01×10−25.01\times 10^{-2}3.35×10−23.35\times 10^{-2}3.29×10−23.29\times 10^{-2}1.71×10−21.71\times 10^{-2}
Table 1: Moment residuals for learned NeMDO and SPH operators. Mean absolute error (MAE) and standard deviation are reported for each monomial moment, averaged over independently perturbed neighbourhood realisations. Residuals quantify the extent to which the learned discrete operators satisfy the imposed Taylor consistency constraints (smaller values are better). The superscripts in the learned operator identify the target operator.
Figure 3: Convergence of the discrete xx -derivative ( left ) and Laplacian ( right ) on a smooth test function with ϵ=0.5\epsilon=0.5 , noting that NeMDO was trained with ϵ=1.0\epsilon=1.0 . Relative L2L_{2} error versus particle spacing ss for the learned NeMDO operator, a second-order LABFM operator, and two uncorrected SPH kernels (quintic spline and Wendland C2).
Figure 4: Normalised eigenvalue spectrum of the global xx -derivative operator with node disturbance ϵ=1.0\epsilon=1.0 and 2500 nodes. Left : xx -derivative operators; right : Laplacian operators.
Figure 5: Resolving power of the x−x- derivative for LABFM, NeMDO, and SPH (with the Quintic Spline and Wendland C2). The left plot shows the real part of the modal response, and the right plot shows the imaginary part on disordered particle distributions with noise level ϵ\epsilon . The horizontal axis shows the true wavenumber kk normalised by the Nyquist wavenumber kNyk_{\mathrm{Ny}} , while the vertical axis shows the corresponding effective wavenumber keff/kNyk_{\mathrm{eff}}/k_{\mathrm{Ny}} . The dashed black line indicates the ideal response, corresponding to a spectral method. The inset on the left plot shows the absolute difference between the different operators and the ideal spectral response. The inset on the right plot shows zoomed in graph of ℑ{k^eff}\Im\{\hat{k}_{eff}\} for range 0≤k^≤0.30\leq\hat{k}\leq 0.3 . The full and dotted lines indicate ky/kx=0k_{y}/k_{x}=0 and ky/kx=1k_{y}/k_{x}=1 , respectively.
Figure 6: Resolving power of the Laplacian operator, the left plot shows the real part of the modal response and the right plot shows the imaginary part on disordered particle distributions with noise level ϵ=1.0\epsilon=1.0 . The horizontal axis shows the true wavenumber qq normalised by the Nyquist wavenumber kNyk_{\mathrm{Ny}} , while the vertical axis shows the corresponding effective wavenumber. The dashed black line indicates the ideal response, corresponding to a spectral method. The inset on the left plot shows the absolute difference between the different operators and the ideal spectral response. The inset on the right plot shows zoomed in graph of ℑ{q^eff2}\Im\{\hat{q}_{eff}^{2}\} for range 0≤q^≤0.30\leq\hat{q}\leq 0.3 . The full and dotted lines indicate ky/kx=0k_{y}/k_{x}=0 and ky/kx=1k_{y}/k_{x}=1 , respectively.
Figure 7: Evaluation of the trade-off between parameter count, wall-clock forward time, and accuracy for NeMDO and baseline kernel methods when computing stencil weights. The vertical axis shows the L2L_{2} error of the xx -derivative of the test function in Equation ( 22 ) at a different resolutions; horizontal lines mark the resolution at which each operator converges (i.e., further refinement does not reduce this error). All NeMDO models are trained with geometric noise ϵ=1.0\epsilon=1.0 , while the convergence tests are performed at ϵ=0.5\epsilon=0.5 .
Figure 8: Velocity magnitude for the Taylor–Green vortex with different operators at different resolutions. Colour axis are rescaled at different time scales for better visualisation of the flow structures.
Figure 9: Relative root mean squared error of different operators with respect to the analytical solution (Equation ( 27 )) for velocity. ss indicates the average particle spacing, and t∗t^{*} indicates characteristic time units of the flow.
∣Ni∣|\mathcal{N}_{i}|OperatorMetricxxyyx2/2x^{2}/2xyxyy2/2y^{2}/2
10NeMDOp=2\text{NeMDO}_{p=2}MAE6.44×10−36.44\times 10^{-3}6.58×10−36.58\times 10^{-3}9.56×10−39.56\times 10^{-3}9.97×10−39.97\times 10^{-3}9.95×10−39.95\times 10^{-3}
St.d.9.89×10−39.89\times 10^{-3}1.05×10−21.05\times 10^{-2}1.76×10−21.76\times 10^{-2}1.70×10−21.70\times 10^{-2}1.63×10−21.63\times 10^{-2}
25NeMDOp=2\text{NeMDO}_{p=2}MAE7.84×10−47.84\times 10^{-4}8.12×10−48.12\times 10^{-4}5.23×10−45.23\times 10^{-4}6.47×10−46.47\times 10^{-4}6.29×10−46.29\times 10^{-4}
St.d.1.06×10−31.06\times 10^{-3}1.10×10−31.10\times 10^{-3}7.21×10−47.21\times 10^{-4}9.13×10−49.13\times 10^{-4}8.91×10−48.91\times 10^{-4}
50NeMDOp=2\text{NeMDO}_{p=2}MAE4.43×10−44.43\times 10^{-4}3.95×10−43.95\times 10^{-4}4.00×10−44.00\times 10^{-4}4.15×10−44.15\times 10^{-4}4.14×10−44.14\times 10^{-4}
St.d.5.68×10−45.68\times 10^{-4}5.18×10−45.18\times 10^{-4}5.25×10−45.25\times 10^{-4}5.76×10−45.76\times 10^{-4}5.60×10−45.60\times 10^{-4}
Table 2: NeMDO moment residuals for xx -derivative operators with varying stencil size ∣Ni∣|\mathcal{N}_{i}| . Mean absolute error (MAE) and standard deviation are reported for each monomial moment, averaged over independently perturbed neighbourhood realisations. Residuals quantify the extent to which the learned discrete operators satisfy the imposed Taylor consistency constraints (smaller values are better).
Figure 10: Parametric results for varying stencil size nn for node disorder of ϵ=1.0\epsilon=1.0 . Colours indicate value of weights. To generate the plots, multiple stencils with normalised relative positions have their weights predicted and plotted on top of each other.
Figure 11: Convergence of the discrete xx -derivative operator on a smooth test function with ϵ=1.0\epsilon=1.0 . Relative L2L_{2} norm versus particle spacing ss for the learned NeMDO operator with varying stencil size. The superscript in the operator indicates the number of particles within the computational stencil.
Figure 12: Resolving power of the gradient operator for varying number of particles in computational stencil. The left plot shows the real part of the modal response, and the right plot shows the imaginary part on disordered particle distributions with noise level ϵ=1.0\epsilon=1.0 . The horizontal axis shows the true wavenumber kk normalised by the Nyquist wavenumber kNyk_{\mathrm{Ny}} , while the vertical axis shows the corresponding effective wavenumber keff/kNyk_{\mathrm{eff}}/k_{\mathrm{Ny}} . The dashed black line indicates the ideal response, corresponding to a spectral method. The inset on the left plot shows the absolute difference between the different operators and the ideal spectral response. The inset on the right plot shows zoomed in graph of ℑ{k^eff}\Im\{\hat{k}_{eff}\} for range 0≤k^≤0.30\leq\hat{k}\leq 0.3 . The full and dotted lines indicate ky/kx=0k_{y}/k_{x}=0 and ky/kx=1k_{y}/k_{x}=1 , respectively.
Figure 13: Normalised third-order NeMDO for the x−x- derivative, colours indicate weight value. Model was trained with ∣Ni∣=50|\mathcal{N}_{i}|=50 and particle disturbance of ϵ=1.0\epsilon=1.0 (and inferred with the same noise).
Metricxxyyx2/2x^{2}/2xyxyy2/2y^{2}/2x3/6x^{3}/6x2y/2x^{2}y/2xy2/2xy^{2}/2y3/6y^{3}/6 
MAE2.08×10−32.08\times 10^{-3}2.07×10−32.07\times 10^{-3}1.71×10−31.71\times 10^{-3}1.92×10−31.92\times 10^{-3}1.85×10−31.85\times 10^{-3}2.42×10−32.42\times 10^{-3}2.43×10−32.43\times 10^{-3}2.39×10−32.39\times 10^{-3}2.86×10−32.86\times 10^{-3} 
St.d.2.70×10−32.70\times 10^{-3}2.73×10−32.73\times 10^{-3}2.24×10−32.24\times 10^{-3}2.52×10−32.52\times 10^{-3}2.40×10−32.40\times 10^{-3}3.29×10−33.29\times 10^{-3}3.21×10−33.21\times 10^{-3}3.21×10−33.21\times 10^{-3}3.61×10−33.61\times 10^{-3} 
Table 3: NeMDO moment residuals for xx -derivative operators with Taylor monomials up to third order, NeMDOp=3x\text{NeMDO}_{p=3}^{x} . Mean absolute error (MAE) and standard deviation are reported for each monomial moment, averaged over independently perturbed neighbourhood realisations. Residuals quantify the extent to which the learned discrete operators satisfy the imposed Taylor consistency constraints (smaller values are better).
Figure 14: Convergence of the discrete xx -derivative operator on a smooth test function with ϵ=0.5\epsilon=0.5 . Relative L2L_{2} norm versus particle spacing ss for the learned NeMDO operator with varying Taylor truncation. The superscript in the NeMDO operators indicate the stencil size ∣Ni∣|\mathcal{N}_{i}| .
Figure 15: Resolving power of the gradient operator for increasing Taylor truncation included during training. The left plot shows the real part of the modal response, and the right plot shows the imaginary part on disordered particle distributions with noise level ϵ=1.0\epsilon=1.0 . The horizontal axis shows the true wavenumber kk normalised by the Nyquist wavenumber kNyk_{\mathrm{Ny}} , while the vertical axis shows the corresponding effective wavenumber keff/kNyk_{\mathrm{eff}}/k_{\mathrm{Ny}} . The dashed black line indicates the ideal response, corresponding to a spectral method. The inset on the left plot shows the absolute difference between the different operators and the ideal spectral response. The inset on the right plot shows zoomed in graph of ℑ{k^eff}\Im\{\hat{k}_{eff}\} for range 0≤k^≤0.30\leq\hat{k}\leq 0.3 . The full and dotted lines indicate ky/kx=0k_{y}/k_{x}=0 and ky/kx=1k_{y}/k_{x}=1 , respectively.
Trained ϵ\epsilonOperatorInferred ϵ\epsilonMetricxxyyx2/2x^{2}/2xyxyy2/2y^{2}/2
1.0NeMDOp=2\text{NeMDO}_{p=2}1.0MAE7.88×10−47.88\times 10^{-4}7.73×10−47.73\times 10^{-4}8.51×10−48.51\times 10^{-4}8.86×10−48.86\times 10^{-4}8.91×10−48.91\times 10^{-4}
St.d.1.07×10−31.07\times 10^{-3}9.89×10−49.89\times 10^{-4}1.21×10−31.21\times 10^{-3}1.24×10−31.24\times 10^{-3}1.22×10−31.22\times 10^{-3}
0.5MAE4.38×10−44.38\times 10^{-4}4.88×10−44.88\times 10^{-4}4.23×10−44.23\times 10^{-4}3.93×10−43.93\times 10^{-4}4.81×10−44.81\times 10^{-4}
St.d.4.74×10−44.74\times 10^{-4}4.78×10−44.78\times 10^{-4}5.26×10−45.26\times 10^{-4}5.04×10−45.04\times 10^{-4}5.79×10−45.79\times 10^{-4}
0.1MAE3.60×10−43.60\times 10^{-4}5.83×10−45.83\times 10^{-4}2.57×10−42.57\times 10^{-4}3.56×10−43.56\times 10^{-4}3.48×10−43.48\times 10^{-4}
St.d.1.03×10−41.03\times 10^{-4}9.63×10−59.63\times 10^{-5}1.76×10−41.76\times 10^{-4}1.43×10−41.43\times 10^{-4}1.66×10−41.66\times 10^{-4}
Table 4: NeMDO moment residuals for xx -derivative operators with varying particle disorder ϵ\epsilon . Mean absolute error (MAE) and standard deviation are reported for each monomial moment, averaged over independently perturbed neighbourhood realisations. Residuals quantify the extent to which the learned discrete operators satisfy the imposed Taylor consistency constraints (smaller values are better). Trained ϵ\epsilon denotes to the particle disorder during training, and inferred ϵ\epsilon denotes the particle disorder during inference.
Figure 16: Convergence of the discrete xx -derivative operator on a smooth test function (Equation ( 22 )) for LABFM p=2 , SPH (with the quintic spline and Wendland C2 kernels) and NeMDO, with varying particle disorder ϵ\epsilon . Full lines indicate inference with ϵ=1.0\epsilon=1.0 , dashed lines indicate inference with ϵ=0.5\epsilon=0.5 , and dotted lines indicate inference with ϵ=0.1\epsilon=0.1 .
Appendix figures & tables2 assets

Supplementary material from the paper’s appendix.

Appendix

Operator∣Ni∣|\mathcal{N}_{i}|# ParametersFhF_{h}MLP Hidden Layers (Emb/Out/Msg/Upd)Graph Layersϵ\epsilon
NeMDO p=2x{}_{p=2}^{x}352.2M256331.0
NeMDO p=2Δ{}_{p=2}^{\Delta}352.2M256331.0
Appendix
Table 5: NeMDO architectures used in the main experiments.
Operator∣Ni∣|\mathcal{N}_{i}|# ParametersFhF_{h}MLP Hidden Layers (Emb/Out/Msg/Upd)Graph Layersϵ\epsilon
NeMDO p=21{}_{p=2}^{1}1011.1k32121.0
NeMDO p=22{}_{p=2}^{2}1546.3k64121.0
Appendix
Table 6: NeMDO architectures used in the cost–accuracy experiments.

Explore similar work

CardsList
May 13, 2026

Topology-Preserving Neural Operator Learning via Hodge Decomposition

In this paper, we study solution operators of physical field equations on geometric meshes from a function-space perspective. We reveal that Hodge orthogonality fundamentally resolves spectral interference by isolating unlearnable topological degrees of freedom from learnable geometric dynamics, enabling an additive approximation confined to structure-preserving subspaces. Building on Hodge theory and operator splitting, we derive a principled operator-level decomposition. The result is a Hybrid Eulerian-Lagrangian architecture with an algebraic-level inductive bias we call Hodge Spectral Duality (HSD). In our framework, we use discrete differential forms to capture topology-dominated components and an orthogonal auxiliary ambient space to represent complex local dynamics. Our method achieves superior accuracy and efficiency on geometric graphs with enhanced fidelity to physical invariants. Our code is available at https://github.com/ContinuumCoder/Hodge-Spectral-Duality

Dongzhe Zheng, Tao Zhong, Christine Allen-Blanchette

Princeton University

Hodge DecompositionPDE Operator LearningNeural Operators
AskPDF
Aug 30, 2026

Joint Spatiotemporal Spectral Neural Operators for Learning PDEs on Irregular Domains

Learning solution operators for partial differential equations (PDEs) on irregular and geometry-dependent domains remains a central challenge in scientific machine learning. While spectral methods provide strong inductive biases for modeling global interactions, they are typically limited to regular domains, and existing neural approaches often require domain warping, interpolation, or costly geometric embeddings. We introduce the Graph Spectral Neural Operator (GSNO), a neural operator that combines spatial graph spectral decompositions with temporal Fourier transforms through a unified space--time spectral kernel. This formulation enables globally coherent operator learning on non-Cartesian discretizations without domain warping or autoregressive rollouts. By replacing learned geometric embeddings with a graph Laplacian spectral basis, GSNO provides geometry-aware spectral learning with low parameter complexity. Across steady and unsteady PDE benchmarks on irregular and geometry-dependent domains, GSNO achieves strong accuracy with reduced runtime and parameter counts, while demonstrating robust zero-shot generalization across mesh resolutions and geometry families.

Abdolmehdi Behroozi, Chaopeng Shen

Department of Civil and Environmental Engineering, Penn State University, University Park, PA 16802, USA

PDE Surrogate ModelingPDE Operator LearningSpectral MethodsNeural Operators
AskPDF
Sep 28, 2026

Preconditioned Physics-Informed Neural Operator Training

Neural operators are typically trained in a supervised fashion, which requires a dataset to be generated with a classical solver. Training them physics-informed, i.e., purely from the governing equations, removes this large offline cost and allows fresh samples to be drawn at every optimization step, but has so far been limited to simplified problems and trails supervised training in accuracy. The obstacle is the ill-conditioning of physics-informed losses, which differential operators induce and which worsens as the discretization is refined. We therefore propose a preconditioned residual loss function and show mesh-independent conditioning for elliptic problems and greatly improved conditioning for saddle point problems. Realized through geometric and algebraic multigrid, the construction applies to linear and nonlinear equations, steady or time-dependent, on structured and unstructured meshes, is agnostic to the neural operator architecture, and adds no cost at inference. On the Poisson, Allen-Cahn and stationary Stokes equations, the resulting label-free training matches supervised training and is four to twenty-five times more accurate than previous physics-informed operator learning methods.

Shizheng Wen, Siddhartha Mishra, Marius Zeinhofer

ETH Zürich

Physics-Informed Neural OperatorsPDE Operator LearningNeural Network Optimization
AskPDF