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

  1. Topology-Preserving Neural Operator Learning via Hodge Decomposition

    May 13, 2026Dongzhe Zheng, Tao Zhong, Christine Allen-BlanchetteHodge DecompositionPDE Operator Learning

    Princeton University

  2. Joint Spatiotemporal Spectral Neural Operators for Learning PDEs on Irregular Domains

    Aug 30, 2026Abdolmehdi Behroozi, Chaopeng ShenPDE Surrogate ModelingPDE Operator Learning

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

  3. Preconditioned Physics-Informed Neural Operator Training

    Sep 28, 2026Shizheng Wen, Siddhartha Mishra, Marius ZeinhoferPhysics-Informed Neural OperatorsPDE Operator Learning

    ETH Zürich