Learning Spectral-Like Mesh-Free Discretisations

arXiv:2609.02833v1 [physics.comp-ph] 2 Sep 2026

Learning Spectral-Like Mesh-Free Discretisations Lucas Gerken Starepravo, Henry Broadley, Jack King

Steven Lind

School of Engineering, University of Manchester Manchester, UK lucas.gerkenstarepravo@postgrad.manchester.ac.uk

School of Engineering, Cardiff University Cardiff, UK

I. I NTRODUCTION Meshfree and particle numerical methods—such as smoothed particle hydrodynamics (SPH) [1], radial basis functiongenerated finite differences (RBF-FD) [2], and the local anisotropic basis function method (LABFM) [3]—are frameworks that can be used to solve partial differential equations without requiring topological connectivity. In multi-scale problems, resolving both the large and the small scales admitted by the discretisation is essential to obtaining accurate solutions. The resolving power—a measure of the accuracy of an operator as a function of the wavenumber—of traditional mesh-free operators is determined mainly by the kernel, the assembly algorithm used to construct the local linear systems, and the order of polynomial consistency imposed [4]. Ensuring polynomial consistency only constrains the operator in the low-wavenumber limit; no mechanism selects for accuracy at the wavenumbers where finescale content actually resides. Whilst spectral methods exactly represent all modes below the Nyquist limit [5], the wavenumber response of polynomial-reconstruction-based methods deviates from the ideal, with errors increasing with increasing wavenumber. Compact (implicit) meshfree schemes [6], [7] can improve resolving power, but couple the derivative values across the domain, requiring the solution of a global linear system at every evaluation. Consequently, features such as dynamic hpadaptivity and non-trivial boundary conditions are challenging to implement in compact schemes. More recently, multi-kernel approaches to optimising resolving power have been proposed [9], improving over traditional explicit operators. These gains are confined to a narrow band wavenumber, however, and the formulation raises open questions: which kernel combination is optimal? How a framework demonstrated in practice with K ≈ 2 behaves when extended to tens or hundreds of kernels? In this work we introduce Spectral-like Neural Discretisation (SpeND), a data-driven approach for obtaining explicit discretisations with spectral-like resolution on unstructured node distributions. We formulate the computation of discretisation weights as a learning problem: a hard-constrained projection layer guarantees that the network outputs satisfy polynomial consistency exactly, and the network learns, within the resulting constrained manifold, the weights that minimise operator error over a prescribed band-limited function space. Training is self-supervised and physics-agnostic, requiring no reference

solutions, and the learned operator is invariant to resolution and translation. We present modified-wavenumber analyses demonstrating that the learned explicit operators achieve spectral-like resolution on significantly disordered node distributions. II. BACKGROUND A. Numerical Approximation of Differential Operators d Let P := {xa }P a=1 ⊂ Ω ⊂ R be an unstructured point cloud sampling a field ϕ. We associate to each point xa the (N ) local neighbourhood Na := {xb ∈ P : ∥xba ∥2 ≤ da }, where (N ) da is the distance from xa to its N -th nearest neighbour, so that |Na | = N . Our objective is to construct, at each node, a discrete approximation LD of a differential operator D acting on ϕ, taking the form of the weighted sum (1) over Na . The entire approximation is thus specified by the stencil weights D }j∈Na , and constructing LD reduces to choosing them. {wb,a X D LD [ϕ(xa )] = ϕba wb,a . (1) j∈Na

Traditional consistent operators such as corrected SPH [8], D RBF-FD [2] and LABFM [3] obtain wb,a by imposing polynomial consistency: the discrete operator must reproduce the action of D exactly on a set of basis functions. Let {pm }M m=1 span the Taylor monomials of degree ≤ q in d dimensions. X D m = 1, . . . , M, (2) pm (xba ) wb,a = D[pm ] x=0 , b∈Ni

which is a linear system Ca wa = b with Ca ∈ RM ×N depending only on the local node geometry. For N > M this system is underdetermined: its solution set is an affine subspace of dimension N −M , and (2) alone does not determine the stencil. Traditional methods close the system implicitly — through the choice of smoothing kernel, basis preconditioning, or a minimum-norm condition. Consistency is therefore imposed with care, whilst the remaining N − M degrees of freedom are fixed as a by-product of that closure. Nothing in this construction refers to the operator’s response at varying wavenumber, and the resolving power of the resulting stencil is whatever the choice of closure happens to deliver. B. Resolving Power Resolving power [10] — the range of wavenumbers over which the discrete operator reproduces the continuous one —

20th SPHERIC World Conference

Rome, Italy, 1-4 February 2027 form an affine subspace of dimension N − M , and we map the raw network output w̃aD onto its closest point,
D D ŵaD = w̃aD − C+ , (5) a Ca w̃a − b

distinguishes two stencils of identical order at the resolutions actually used, and can be quantified by modal analysis. Since LD is linear, it acts on Fourier modes independently, so substituting eik·xba into (1) with node a at the origin characterises the operator completely. Taking D = ∂x , for which ∂x eik·xba = ikx eik·xba , we define the effective wavenumber keff ≡ −i L∂x [eik·xba ], normalised so that the exact operator recovers keff = kx : X ∂x wb,a sin(k·xba ) + keff =

with C+ a the Moore–Penrose pseudoinverse. Equation (5) is the minimum-norm correction onto the solution set of (2), reducing to the identity whenever w̃aD is already consistent. The projection places no requirement on the architecture that computes w̃aD . IV. R ESULTS We report preliminary results with fθ instantiated as a graph neural network. For each neighbourhood Na we define a graph Ga := (Va , Ea ), with node set Va ≡ Na comprising particle a and its neighbours, and edge set Ea := {(a, b), (b, a) : b ∈ Na \ {a}} giving star-shaped, bidirectional connectivity to the central particle. The network is trained on two-dimensional disordered node distributions, with a hard-constraint layer imposing fourthorder consistency for ∂x (so b = [1, 0, . . . , 0]⊤ ∈ R14 ), and the modal target of (4) set to the exact response cos(k · xba ) over ∥k∥ ∈ (0, kny ).

b∈Na

i

X

∂x wb,a [1 − cos(k · xba )] .

(3)

b∈Na

The deviation of ℜ{keff } from kx is dispersion error, mispropagating under-resolved modes; ℑ{keff } is dissipation error, with positive values implying unconditional instability. III. M ETHOD In SpeND, the discretisation weights are parametrised by a neural network conditioned on the local stencil geometry D {ŵb,a }j∈Na = fθ ({x̂ba }b∈Na ), where θ are the network parameters and are shared across all stencils and trained once, and x̂ba = xba / maxb∈Na ∥xba ∥2 normalises each stencil to the unit D ball, and ŵb,a represents weights normalised by the local stencil size, rendering fθ invariant to the physical scale of the cloud; the scale factor is reintroduced at inference.

A. Modal Analysis Fig. 1 compares the effective wavenumber of the learned operator against explicit LABFM at equal stencil size |Na | = 30 on the same node unstructured distribution, and against a 4thorder finite difference scheme on a structured grid. Two wave angles are shown, ky = 0 (solid) and ky = kx (dotted), with wavenumbers normalised by kny . SpeND follows the spectral response over a substantially wider range of wavenumbers than either LABFM or FDM at equal stencil size, departing only close to the Nyquist limit. The inset shows absolute error over k/kny ∈ (0, 0.6): SpeND outperforms both methods across the full band, not only near the Nyquist limit. These results demonstrate: (i) the benefit of optimising the modal response directly, with SpeND resolving substantially more of the band than either baseline at equal stencil size; (ii) that this is achieved without sacrificing consistency

A. Training The analysis of §II-B makes “spectrally accurate” concrete: for every k in the resolvable band, keff should coincide with the exact response, which for D = ∂x is cos(k · xba ). SpeND trains against this criterion directly, evaluating (3) on the projected weights. Let K tile the half-disk {∥k∥ ≤ kny , kx ≥ 0}; the omitted half follows from the parity of (3). The per-stencil loss penalises dispersion and dissipation error, 1 X 2 2 Li = ω(k) (ℜ{keff } − kx ) + γ (ω+ (k)ℑ{keff }) , (4) |K| with the total objective L(θ) = Ea∼P [Li ] averaged over the stencils of a pre-generated set of disordered point clouds, and γ > 0 balancing resolving power against stability. Both weights are asymmetric: ω = λ> > 1 where ℜ{keff } > kx and unity otherwise, since a phase speed below the exact one merely under-resolves a mode whereas one above it propagates energy ahead of the physical wavefront; likewise ω+ = λ+ > 1 where ℑ{keff } > 0, improving stability of the discretisation weights.

1.0

̂ } ℜ{keff

0.8

|ℜ{k̂ eff} − k|̂

k∈K

10−3 10−6 10−9

0.6

0.2

0.4

SpeND (ours)

0.6

LABFM FDM

0.4

Spectral

0.2 0.0 0.0

B. Hard Constraints Rather than penalising consistency in the loss like in [11], SpeND enforces it exactly, so that every operator the network can express satisfies (2) by construction and the learned degrees of freedom are precisely the N − M directions left free by it. Since N > M and Ca has full row rank, the consistent weights

0.2

0.4

0.6

0.8

1.0

Fig. 1. Effective wavenumber ℜ{keff } of the ∂x operator against k, both normalised by kny , for SpeND and explicit LABFM at p = 4 with |Na | = 30 on disorder node distributions, and a 4th-order finite difference scheme on a structured grid. Solid lines denote ky = 0, dotted ky = kx ; the diagonal is the exact response. Inset: absolute error |ℜ{keff } − kx | over k/kny ∈ (0, 0.6).

2

20th SPHERIC World Conference

Rome, Italy, 1-4 February 2027

— the hard-constrained layer leaves the fourth-order asymptotic behaviour intact, visible as the rapid decay of absolute error as k → 0, so the resolving power is not bought at the cost of formal order; (iii) the effectiveness of the asymmetric dispersion penalty, with no overshoot of the exact response anywhere in the band.

The operator is explicit, and so requires no global solve at inference, unlike compact schemes. In Eulerian simulations the weights are computed once and reused, making SpeND identical in cost to SPH, LABFM or RBF-FD at run time; thus, the resolving power is obtained for no additional expense. This is particularly relevant for non-linear problems such as the Navier– Stokes equations, where the upper third of the spectrum must be dealiased regardless. Fig. 1 shows SpeND to be near-exact across the entire range of wavenumbers that survive dealiasing, whereas both baselines depart from the exact response well inside it. Work in progress addresses the remaining validation: out-ofdistribution validation, and time-dependent solves demonstrating eigenvalue stability and sustained resolving power over long integration times.

B. Convergence Fig. 2 shows the L2 error of the ∂x operator under refinement of the node spacing s, on unstructured node distributions, for the test field ϕ(x, y) =

7 X

4 1 sin(2πy) sin 2πj x − 41 , π j j=1

(6)

j odd

ACKNOWLEDGMENT J. K. is funded by the Royal Society via a University Research Fellowship (URF\R1\221290). We would like to acknowledge the assistance given by Research IT and the use of the Computational Shared Facility at the University of Manchester.

a four-term truncated square wave with harmonics up to k = 14π. Both schemes recover the expected fourth-order rate once the highest harmonic is resolved, as guaranteed for SpeND by the hard-constrained layer of §III-B. SpeND is more accurate throughout, by a factor of approximately 2 at the coarsest resolution widening to 3 at the finest. SpeND recovers the classical convergence behaviour of a p = 4 discretisation whilst resolving substantially more of the band at practical resolutions. Note also that (6) is a superposition of harmonics with non-uniform amplitudes rather than a single mode of the form appearing in the training objective.

R EFERENCES [1] L. B. Lucy, “A Numerical Approach to the Testing of the Fission Hypothesis,” Astron. J., vol. 82, pp. 1013–1024, Dec. 1977. [2] V. Shankar, G. B. Wright, R. M. Kirby, and A. L. Fogelson, “A Radial Basis Function (RBF)-Finite Difference (FD) Method for Diffusion and Reaction-Diffusion Equations on Surfaces,” arXiv preprint arXiv:1404.0812, April 2014. [3] J. R. C. King, S. J. Lind, and A. M. A. Nasar, “High Order Difference Schemes Using the Local Anisotropic Basis Function Method,” J. Comput. Phys., vol. 415, p. 109549, Aug. 2020. [4] J. R. C. King and S. J. Lind, “High-Order Simulations of Isothermal Flows Using the Local Anisotropic Basis Function Method (LABFM),” J. Comput. Phys., vol. 449, p. 110760, Jan. 2022. [5] A. Brandenburg, “Computational aspects of astrophysical MHD and turbulence,” in Advances in Nonlinear Dynamos, CRC Press, 2003, pp. 269–344. [6] N. Trask, M. Maxey, and X. Hu, “Compact moving least squares: An optimization framework for generating high-order compact meshless discretizations,” J. Comput. Phys., vol. 326, pp. 596–611, Dec. 2016. [7] H. M. Broadley, S. J. Lind, and J. R. C. King, “Compact LABFM: a framework for meshless methods with spectral-like resolving power,” arXiv preprint arXiv:2603.11668, March 2026. [8] J. Bonet and T.-S. L. Lok, “Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations,” Comput. Methods Appl. Mech. Eng., vol. 180, no. 1, pp. 97–115, Dec. 1999. [9] H. Broadley, J. R. C. King, and S. J. Lind, “Improving the accuracy of meshless methods via resolving power optimisation using multiple kernels,” arXiv preprint arXiv:2510.20365, Oct. 2025. [10] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. Comput. Phys., vol. 103, no. 1, pp. 16–42, Nov. 1992. [11] L. G. Starepravo, G. Fourtakas, S. Lind, A. B. Harish, T. Tang, and J. R. C. King, “Learning mesh-free discrete differential operators with self-supervised graph neural networks,” arXiv preprint arXiv:2603.24641, March 2026.

V. C ONCLUSIONS We have presented SpeND, a framework for learning explicit mesh-free discretisation operators with spectral-like resolving power on unstructured node distributions. By enforcing polynomial consistency exactly through a hard-constrained projection layer, and optimising the modal response of the operator within the remaining degrees of freedom, SpeND resolves substantially more of the wavenumber band than explicit LABFM and fourthorder finite differences at equal stencil size, whilst retaining the fourth-order convergence rate of the imposed consistency.

L2 norm

10−1

10−2

LABFM SpeND(ours)

10−3

10−4

10−2

s

Fig. 2. Convergence of the ∂x operator under refinement of the node spacing s, for SpeND and explicit LABFM at p = 4 with |Na | = 30, on unstructured node distributions. Error is the L2 norm over all nodes of the field (6). Both schemes attain fourth order; SpeND is a factor of 2–3 more accurate throughout.

3

Related documents

• Learning Spectral-Like Mesh-Free Discretisations#924956
• Learning Spectral-Like Mesh-Free Discretisations#920722
• ADELIA: Automatic Differentiation for Efficient Laplace Inference Approximations#168263
• PTL-Diffusion: Manifold-Aware Diffusion with Periodic Terminal Laws#267679
• Weighted Laplacian Flow: A Deterministic Particle Flow with Exponential Convergence#1084330