Learnable Splines for Unstructured Neural Operators
Authors: Vignesh Gopakumar, Zhuoyuan Wang, Ander Gray, Valentin Duruisseaux, Josh Williams, William Michael Edward Ellis, Nuno Nobre, Matthew Field, James Panayis, Alexander I. Blair, Andrew Davis, Stanislas Pamela, Robert Akers, Marc Peter Deisenroth, Anima Anandkumar
Organizations: Centre for Artificial Intelligence, University College London · Computing + Mathematical Sciences, Caltech · LIX, École Polytechnique, CNRS, Institut Polytechnique de Paris · Hartree Centre, Science and Technology Facilities Council · Computing Division, UK Atomic Energy Authority
Neural operators approximate solution operators of partial differential equations at a fraction of the cost of numerical solvers. Efficient grid-based architectures such as the Fourier neural operator (FNO) assume uniform rectangular grids, whereas engineering simulations use large, irregular, and locally refined meshes. We introduce L-SplineFNO, a memory-efficient neural operator that couples an FNO backbone with a sparse, learnable tensor-product B-spline transfer between unstructured meshes and a latent Cartesian grid. A local normalised scatter avoids the geometry-dependent solve of analytic spline transformation while preserving compact support, partition of unity, and spline smoothness by construction. The neural operator, conditioned on global design and operating parameters, learns the mapping of geometry to the field solution. Across six three-dimensional engineering benchmarks, L-SplineFNO achieves strong predictive accuracy while requiring substantially less memory and compute than baseline operators, scaling to 26.6 million nodes in training and 46 million in inference on a single A100 GPU. It also supports zero-shot super-resolution and efficient closed-form spatial derivatives. Uniquely, the analytic spline decoder yields rigorous prediction bounds over continuous space and uncertain parameters, via Bernstein convex hulls and Taylor models, demonstrated on a 29.2-million-node MAST-U tokamak model.
Figures & tables
Figure 1: End-to-end conditional L-SplineFNO architecture. The unstructured geometry and global conditioning vector are first mapped to conditioned nodal features through coordinate normalisation, positional encoding, and FiLM modulation. A learnable Cp−1 spline scatter then transforms these features from the unstructured mesh to a structured latent grid, where they are processed by a Fourier neural operator (FNO) backbone. The resulting grid representation is mapped back to the physical mesh using analytic spline reconstruction and a pointwise projection to produce the predicted field. This construction provides a sparse, structure-preserving, differentiable mesh-to-grid-to-mesh mapping with the structured Fourier processing.
Figure 2: Error maps produced by baseline architectures in modelling heat distribution across the heat sink. L-SplineFNO produces the most accurate predictions, with minimal errors concentrated at the edges of the heat sink where the temperature gradient is highest.
Model
Heat Sink
Bracket
Aircraft
Building Fire
Heated Duct
HIVE Monoblock
N=16,384
N=8,192
N=131,072
N=98,496
N=124,800
N=194,691
Baselines
Transolver
0.040±0.002
0.591±0.025
8.32±0.88
7.341±0.220
0.473±0.050
0.301±0.061
00:41
01:48
01:26
01:17
01:15
02:28
UPT
0.039±0.001
0.633±0.014
13.87±0.43
7.162±0.158
0.523±0.053
0.299±0.109
00:33
01:30
00:58
01:00
0:50
01:25
Table 1: Relative L2 error (%, mean ± std over five seeds; top line of each entry) and training time (hh:mm on a single NVIDIA A100; grey line beneath) for each architecture under FiLM conditioning. The best mean error for each case is shown in bold. N is the number of spatial training nodes. Heat Sink, Heated Duct and HIVE predict absolute temperature, so their relative errors are small by construction. L-SplineFNO outperforms baselines on most experiments, while remaining significantly compute-efficient.
J(θ)=M1∑i=1My(i)2y(i)−y(i)2
(12)
Table 4
Figure 3: Error maps produced by baselines in modelling heat distribution across the HIVE monoblock. L-SplineFNO models the temperature across the structure with 99.9% accuracy (sample-specific error shown).
Figure 4: Derivative-evaluation cost and scalability. ( 4(a) ) Mean derivative runtime (left) and additional peak GPU allocation (right), averaged across the benchmark experiments. Forward preparation is excluded; memory is measured above the prepared live state, not total GPU use. Full-model autograd VJP measures coordinate sensitivity, whereas decoder autograd and analytic spline differentiation evaluate pointwise field gradients. Closed-form spline derivatives give exact gradients at a fraction of the cost. ( 4(b) ) Mean largest successful node count across three seeds on an A100 80GB GPU in FP32 with batch size one. Training includes a complete Adam step. L-SplineFNO accomodates 26.6M nodes in training and 46.0M in inference, demonstrating superior scalability against the baselines.
Figure 5: MAST-U virtual tokamak: Displacement magnitude (mm) from the MOOSE simulation (ground truth) and the L-SplineFNO (prediction) and the absolute error across them (mm).
Figure 6: Rigorous enclosure of L-SplineFNO over the uncertain parameter box B , obtained combining a Taylor-remainder bound with the convex-hull property of the spline.
Appendix figures & tables39 assets
Supplementary material from the paper’s appendix.
Appendix
Symbol
Meaning
Problem and data
Ω⊂Rd , d
Spatial domain and its dimension ( d=3 throughout).
x∈Ω
Generic continuous point; written x in the enclosure formulas.
μ∈P⊂Rk
Conditioning vector of design, load, material and operating parameters.
U , u0
Solution space ( Equation 1 ); optional per-node input field, absent in every benchmark ( cin=0 ).
correction index in κθ(ℓ) ; displacement component (enclosures)
k
dimμ
Bernstein multi-index (enclosures)
q
query point qj
collocation index uq ( Section B.2 )
p
spline degree
static pressure (aircraft)
c
model width
spline coefficient ci(μ) (enclosures)
Appendix
Table 3: Symbols with more than one meaning. The context, given in the last column, always fixes which is meant.
Depth n
Partitions 2n
Median width (mm)
Max width (mm)
Time (s)
0
1
8.70
148.12
5.71
1
2
7.65
142.32
10.64
2
4
6.62
130.71
20.34
3
8
4.24
93.26
39.61
4
16
3.89
79.88
81.99
6
64
2.71
57.98
322.43
Appendix
Table 4: MAST-U spatial subdivision ablation at nominal parameters. Depth zero uses no subdivision within each original spline cell; each subsequent depth doubles the number of sub-boxes per cell. Widths are displacement-magnitude enclosure widths over the same 4,077,312 displayed triangles, with equal weight per triangle. Times are medians of three warm, synchronised bound calculations on one NVIDIA A100 GPU, excluding model loading, source encoding and reporting. Depth three is the current subdivision setting.
Figure 7: Displacement-magnitude enclosure widths at fixed face centres for eval_02 over the 1× parameter box, using ordinary FP32, ordinary FP64 and validated binary64. The FP32/FP64 calculations use verified norms and quadratic extrema with ordinary propagation; validated binary64 encloses the real stored-weight surrogate. All panels use the manuscript camera, the same geometry and a shared linear colour scale without clipping. The distributions correspond to Table 5 ; panel annotations give the median and maximum widths over all 4,077,312 face centres.
Interval implementation
Time (s)
Median width (mm)
Max width (mm)
Original formulation, FP32
74.63
2.43694
6.34944
Verified norms/extrema, FP32
106.68
2.43694
6.34944
Verified norms/extrema, FP64
133.83
2.43694
6.34944
Validated binary64
356.40
2.47411
6.45884
Appendix
Table 5: Parameter-box arithmetic comparison for the evaluation case at 1× , using the optimized validated binary64 implementation. Times are median per-box costs over three repetitions in the same A100 allocation with cached constants, excluding setup and warmup. Widths are displacement-magnitude enclosure widths at fixed face centres. The ordinary FP32/FP64 rows have no end-to-end rounding certificate; validated binary64 encloses the real surrogate defined by the stored weights.
Figure 8: FP32 displacement-magnitude enclosure widths for eval_02 over the 0.5× , 1× and 2× parameter boxes in Table 6 . Both parameter half-widths scale together about the same nominal parameters, with no parameter subdivision. The geometry and camera are identical to Figure 7 . A shared linear colour scale spanning the largest width in these three panels shows the growth of field variation and enclosure conservatism without clipping or independent panel rescaling. Widths are evaluated at fixed face centres using the updated ordinary FP32 formulation.
Box factor
α half-width (K−1)
Tw half-width (K)
Median width (mm)
Max width (mm)
0.5×
5×10−8
2.5
0.47310
1.02662
1×
10−7
5
2.43694
6.34944
2×
2×10−7
10
20.55161
51.67701
Appendix
Table 6: FP32 displacement-magnitude enclosure widths for eval_02 across three parameter-box sizes. Both half-widths scale by the box factor about the same nominal parameters. Median and maximum widths are computed over all 4,077,312 fixed face centres. These are ordinary-arithmetic enclosure calculations, without an end-to-end rounding certificate.
Interval implementation
Constant setup/check (s)
Warmup (s)
Original formulation, FP32
2.64
88.75
Verified norms/extrema, FP32
336.39
108.64
Verified norms/extrema, FP64
0.46
134.67
Validated binary64
0.41
522.53
Appendix
Table 7: One-time constant preparation/checks and one full warmup per method, excluding the separately reported restore/geometry cost. The updated methods ran in order FP32, FP64, validated binary64 in one process: FP32 builds the shared verified constants, while the last two rows report cache checks. Warmup uses the 0.5× box; the measured comparison in Table 5 uses 1× .
Dataset
Physics / target
Cases train/test
Native nodes
Used nodes
dimμ
Mesh structure
Heat sink
Conduction / T
625 ( 438/187 )
4,671 – 32,645
16,384
4
Sample-dependent
Bracket
Elasto-plasticity / σ
3,000 ( 2,100/900 )
5,054 – 15,210
8,192
4
Sample-dependent
Aircraft
Compressible flow / (ρ,u,v,w,p)
150 ( 105/45 )
331,971
131,072
3
Shared correspondence
Building fire
Fire LES (time-avg.) / (ΔT,ρair,ρprod)
547 ( 383/164 )
116,928
98,496
7
Sample-dependent, three-storey
Heated duct
Conjugate heat transfer / T
122 ( 103/19 )
124,800
124,800
5
Shared, three-region
HIVE
Induction heating / T
186 ( 149/37 )
194,691
194,691
6
Shared
Appendix
Table 8: Overview of the six benchmarks. “Native nodes” denotes the original mesh size, whereas “used nodes” is the spatial budget presented to the model. Variable-size meshes are sampled or duplicate-padded to this budget.
Figure 9: Distribution of the target field over all nodes of all cases, one panel per benchmark (pressure shown for aircraft and temperature change for building fire). Counts are normalised to a fraction of nodes and plotted on a log ordinate; the horizontal axis is trimmed to [p0.05,p99.95] and the lighter shading marks the tails outside [p1,p99] . The dashed line is the median. Sample skewness (γ1) and kurtosis (β2) are annotated (a Gaussian has Skew=0 , Kurt=3 ).
Parameter
Symbol
Levels
Meaning
Sink width
W
20,22.5,25,27.5,30 mm
y extent; fins span it fully
Fin height
H
15,16.25,17.5,18.75,20 mm
z extent above the baseplate
Thickness
t
1.5,1.75,2,2.25,2.5 mm
fins and baseplate; z spans [−t,H]
Fin count
nfin
4,5,6,7,8
plates stacked along the 26 mm base
Appendix
Table 9: Heat sink: the 54 factorial design. All 625 combinations are present exactly once.
Figure 10: Error maps produced by various architectures in modelling heat distribution across the heat sink.
Figure 11: Error maps produced by various architectures in modelling stress across the bracket.
Channel
min
max
mean
std
Skew
Kurt
per-condition mean (%)
ρ
−0.017
18.24
1.315
0.867
2.33
16.0
35.95
u
−0.324
7.670
5.665
1.957
−1.42
3.20
5.25
v
−3.044
3.500
0.245
0.522
−0.13
4.30
67.77
w
−3.900
4.070
0.063
0.714
−0.05
5.00
61.33
p
−0.035
46.46
1.650
2.131
8.17
107.8
46.73
joint
−3.900
46.46
1.788
2.471
2.08
14.3
23.30
Appendix
Table 10: Aircraft: per-channel statistics over all 150 cases at the 131,072 sampled nodes. All quantities are freestream-normalised and dimensionless. The last column is the relative L2 error of predicting the mean field of the matching flight condition, i.e. the score of a predictor that uses μ perfectly and ignores the geometry; a model must beat 23.3% on the joint field before it has learned anything from the airframe shape.
Figure 12: Error maps produced by various architectures in modelling the pressure distribution across the aircraft.
Parameter
Symbol
Range
Discretisation
Group
Visible in the mesh
Burner centre, x
xf
[2.25,9.75] m
0.25 m
fire
yes, exactly
Burner centre, y
yf
[2.25,9.75] m
0.25 m
fire
yes, exactly
Ventilation rate
Q
[−4.0,4.0] m 3 s -1
continuous
ventilation
no
Left vent centre, y
yL
[1.0,11.0] m
0.25 m
ventilation
no
Left vent centre, z
zL
[4.75,6.5] m
0.25 m
ventilation
no
Right vent centre, y
yR
[1.0,11.0] m
0.25 m
ventilation
no
Appendix
Table 11: Building fire: the seven conditioning parameters, their sampling ranges, and positional discretisation for the 547-case corpus. Left and right denote the vents at x=0 and x=12 m, respectively. The last column indicates whether each parameter is recoverable from the node coordinates; see the redundancy discussion below.
Cell dimensions [m]
Cells
Volume [m 3 ]
Where
0.1253
18,432 (15.8%)
36 (2.5%)
refinement box around the burner
0.25×0.25×0.125
17,280 (14.8%)
135 (9.4%)
z=3 – 4 m excluding refinement box
0.253
81,216 (69.5%)
1269 (88.1%)
everywhere else
Appendix
Table 12: Building fire: the three cell classes of the 116,928-cell mesh (1440 m 3 in total, including solid obstructions), before solid-cell masking and refinement-box removal. On this full Cartesian mesh, an unweighted node average over-weights the refinement box by ∼6.3× relative to its volume, so integral quantities of interest over the gas domain must exclude solid cells and be volume-weighted; the loader carries the cell volume per node for this purpose.
Channel
min
max
mean
std
ΔT [K]
0
189.4
39.32
16.62
ρair [kg m -3 ]
0.763
1.630
1.247
0.0845
ρprod [kg m -3 ]
0
0.1007
0.0462
0.0164
Appendix
Table 13: Building fire: per-channel statistics over all 547 cases at the 92,837 retained gas nodes (solid cells and refinement box removed).
Figure 13: Error maps produced by various architectures in modelling temperature change across the building due to a fire in the lower floor.
Parameter
Symbol
Physical range
Transform
Shape
Heater power density
qh
2.5×106 – 2.0×107 W m -3
log(qh/107)
no
Inlet velocity
Uin
5×10−4 – 2×10−3 m s -1
Uin/10−3
yes
Inlet temperature
Tin
280 – 320 K
(Tin−300)/78.5
no
Heater conductivity
kh
40 – 120 W m -1 K -1
kh/80
yes
Metal conductivity
km
100 – 300 W m -1 K -1
km/200
yes
Appendix
Table 14: Heated duct: the five conditioning parameters. μ is the dimensionless transform, not the raw value — only the log on qh , which spans a factor of eight, is non-trivial. The last column records whether the parameter changes the shape of the temperature field or only rescales it; see the degeneracy noted below.
Figure 14: Error maps produced by various architectures in modelling the temperature distribution from the conjugate heat transfer across the heated duct.
Parameter
Group
Prior
Realised range
Role
nusselt_scale
HTC
U(0.001,0.031)
[0.00112,0.03090]
C in Equation 50
nusselt_re_power
HTC
U(0.8,0.9)
[0.80043,0.89950]
a in Equation 50
nusselt_pr_power
HTC
U(0.35,0.45)
[0.35048,0.44997]
b in Equation 50
coil_target_offset_x
geometry
U(−0.01,0.01) m
[−0.00996,0.00995]
coil position
coil_target_offset_y
geometry
U(−0.01,0.01) m
[−0.00994,0.00996]
coil position
coil_target_offset_z
geometry
U(0.019,0.039) m
[0.01901,0.03892]
coil standoff
Appendix
Table 15: HIVE: the six uncertain parameters, their sampling priors and the ranges realised over the 186 completed runs.
Figure 15: Error maps produced by various architectures in modelling the temperature distribution from inductive heating in a monoblock in the HIVE case.
Parameter
Symbol
Unit
Realised range
Mean
Role
Thermal expansion coefficient
α
10−5 K -1
1.015 – 1.982
1.50
eigenstrain scale
Wall temperature
Tw
K
508.6 – 996.2
750.1
thermal boundary condition
Appendix
Table 16: MAST-U: the two conditioning parameters over the pooled 70 solves. Both are min–max normalised per variable using the pooled corpus ranges before entering the conditioning module.
Model
Encoder
Processor
Decoder
Attention Cost
Transolver
feature slicing
softmax MHSA
deslicing
O(Nnsc+ns2c)
GNOT
trunk/branch MLPs
linear cross + self-attn.
MLP
O(Nc2)
UPT
Perceiver pooling
softmax MHSA
Perceiver query
O(Nnsc+ns2c)
MSPT
patch pooling
chunked MHSA + supernodes
MLP
O(N2/np+Nnpnq)c
NUNO
K-D tree + scatter
FNO on subdomain grid
trig. interpolation
—
GINO
radius-graph GNO
FNO on latent grid
radius-graph GNO
—
Appendix
Table 17: The six baselines by encode–process–decode role. N is the number of mesh nodes, ns the number of slices or latent tokens, np the number of patches and nq the supernodes per patch, c the model width. The two gridded operators have no attention term; their cost is dominated by the FFT on the latent grid, O(GlogG) for G grid points.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,433,413
0:41
0.0589
0.0460
1.2804
0.9989
Transolver
Concatenation
2,434,437
0:41
0.0458
0.0341
1.3449
0.9993
Transolver
Cross-attention
2,464,422
0:44
0.0438
0.0330
1.3281
0.9994
Transolver
FiLM
2,450,503
0:41
0.0385
0.0280
1.3742
0.9995
UPT
None
2,252,097
0:33
0.0907
0.0741
1.2230
0.9972
UPT
Concatenation
2,252,993
0:32
0.0423
0.0308
1.3749
0.9994
Appendix
Table 18: Performance metrics of various models and conditioning methods for the Heat Sink (3D, steady-state, thermal) case. Best value per metric within each architecture in bold.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,433,413
1:47
23.6733
22.9590
1.0311
0.6963
Transolver
Concatenation
2,434,437
1:47
0.9006
0.6526
1.3799
0.9994
Transolver
Cross-attention
2,464,422
1:57
0.8573
0.6495
1.3198
0.9995
Transolver
FiLM
2,450,503
1:49
0.5698
0.4059
1.4039
0.9998
UPT
None
2,252,097
1:28
27.0258
25.6378
1.0541
0.6351
UPT
Concatenation
2,252,993
1:28
1.9136
1.3180
1.4519
0.9975
Appendix
Table 19: Performance metrics of various models and conditioning methods for the Bracket (3D, steady-state, elasto-plastic) case. Best value per metric within each architecture in bold.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,434,441
1:26
64.9073
54.6929
1.1868
0.3622
Transolver
Concatenation
2,435,209
1:27
7.6032
5.7802
1.3154
0.9514
Transolver
Cross-attention
2,464,874
1:36
14.8217
11.5540
1.2828
0.8201
Transolver
FiLM
2,451,467
1:26
9.2876
6.9728
1.3320
0.9273
UPT
None
2,252,997
0:57
62.9624
52.7489
1.1936
0.3524
UPT
Concatenation
2,253,669
0:56
14.9982
11.6138
1.2914
0.8207
Appendix
Table 20: Performance metrics of various models and conditioning methods for the Aircraft (3D, steady-state, compressible flow, 5 fields) case. Best value per metric within each architecture in bold. R2 is the unweighted mean of the per-variable field R2 over the five outputs.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,433,927
1:11
14.6908
11.7212
1.2534
0.4330
Transolver
Concatenation
2,435,719
1:11
7.7780
5.6677
1.3723
0.9355
Transolver
Cross-attention
2,467,432
1:17
8.4583
6.2634
1.3504
0.9195
Transolver
FiLM
2,451,209
1:12
7.4685
5.5415
1.3477
0.9378
UPT
None
2,252,547
0:54
16.7382
13.1177
1.2760
0.4039
UPT
Concatenation
2,254,115
0:54
8.6675
6.2352
1.3901
0.9217
Appendix
Table 21: Performance metrics of various models and conditioning methods for the FDS Building Fire (3D, time-averaged LES, three output fields) case. Single run per entry (seed 42); training times on a single NVIDIA H100; N=92,837 nodes. Best value per metric within each architecture in bold.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,433,413
1:14
18.8793
18.2311
1.0356
0.0207
Transolver
Concatenation
2,434,693
1:13
0.7063
0.5926
1.1918
0.9982
Transolver
Cross-attention
2,465,126
1:22
0.9252
0.8388
1.1029
0.9954
Transolver
FiLM
2,450,567
1:14
0.4481
0.3704
1.2097
0.9993
UPT
None
2,252,097
0:48
18.8869
18.2296
1.0361
0.0181
UPT
Concatenation
2,253,217
0:48
1.1621
0.9795
1.1864
0.9959
Appendix
Table 22: Performance metrics of various models and conditioning methods for the Heated Duct (3D, steady-state, multi-region conjugate heat transfer) case. Best value per metric within each architecture in bold.
Model
Cond. Method
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Transolver
None
2,119,957
2:23
6.3799
5.5841
1.1425
0.2876
Transolver
Concatenation
2,121,397
2:24
0.5016
0.4167
1.2038
0.9953
Transolver
Cross-attention
2,219,766
2:42
0.3052
0.2462
1.2396
0.9983
Transolver
FiLM
2,178,903
2:26
0.2884
0.2463
1.1708
0.9984
UPT
None
2,382,641
1:15
6.3255
5.5387
1.1420
0.2766
UPT
Concatenation
2,384,033
1:16
0.9408
0.8203
1.1470
0.9844
Appendix
Table 23: Performance metrics of various models and conditioning methods for the HIVE Monoblock (3D, steady-state, thermal) case. Best value per metric within each architecture in bold.
Figure 16: Runtime (top) and peak allocated GPU memory (bottom) as node count increases, separated by execution phase. Both axes are logarithmic. Points summarise successful seeds only, and lines connect measured sizes without extrapolation. NUNO’s training result at ten million nodes includes two successful seeds; the remaining seed runs out of memory. Missing points indicate unsuccessful scaling rather than zero cost.
Figure 17: Runtime and peak allocated memory per million nodes, using the same measurements and aggregation as Figure 16 . Both axes are logarithmic. A nearly horizontal curve indicates approximately linear growth of the total cost over the measured range; a decreasing curve can reflect amortisation of fixed overhead. L-SplineFNO maintains a broadly stable cost per node as the point cloud grows.
Experiment
Nodes
Rel. L2 (%)
R2
Heat Sink
31,313
0.0239
0.9996
Bracket
14,703–15,210
3.5142
0.9958
Aircraft
331,971
8.3237
0.9899
Appendix
Table 24: Zero-shot super-resolution with FiLM-conditioned L-SplineFNO.
Figure 18: Accuracy–cost trade-offs for multiscale scattering and spline degree. Values are percentage changes relative to the single-scale and cubic configurations, respectively, computed against the baseline trained at the same seed. Points are means over three seeds and error bars one standard deviation; spline orders 6 and 8 were run at a single seed and carry no bar. Spline orders 2 , 4 , 6 , and 8 correspond to degrees p=3 , 5 , 7 , and 9 , i.e. a Cp−1 basis. By increasing spline order and the scatter length, the model only makes marginal gains against a linear climb of computational costs. Absolute results are reported in Tables 25 and 26 .
Case
Scatter Scales
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Heat-sink
1
2,395,775
0:43
0.02842 ± 0.00331
0.02201 ± 0.00269
1.291
0.99968
1–2
2,396,969
0:51
0.03012 ± 0.00311
0.02303 ± 0.00257
1.309
0.99966
1–2–4
2,397,563
0:59
0.03036 ± 0.00327
0.02310 ± 0.00243
1.314
0.99965
Bracket
1
2,387,731
0:54
1.9825 ± 0.0207
1.1096 ± 0.0283
1.787
0.99832
1–2
2,388,344
1:16
1.9734 ± 0.0284
1.0996 ± 0.0316
1.795
0.99833
1–2–4
2,388,651
1:37
2.0134 ± 0.0328
1.1319 ± 0.0225
1.779
0.99826
Appendix
Table 25: Multiscale L-SplineFNO (FiLM conditioning) across the benchmark cases at increasing scatter-scale depth, all at spline degree 3. Heated-duct and FDS run under region mode basis . Errors, L2/L1 and R2 are mean ± s.d. over seeds 42, 43 and 44 where a cell has several (train/test split and initialisation vary with the seed); training time is the seed mean. For the multi-variable cases (Aircraft, FDS) R2 is the mean of the per-variable scores.
Case
Spline Degree
# Parameters
Train Time (hrs:mins)
Rel. L2 % ( ↓ )
Rel. L1 % ( ↓ )
L2/L1
R2 ( ↑ )
Heat-sink
3
2,395,775
0:42
0.02867 ± 0.00184
0.02215 ± 0.00157
1.295
0.99968
5
2,395,769
1:10
0.03425 ± 0.00750
0.02627 ± 0.00551
1.302
0.99957
7
2,395,778
2:00
0.02852 (+0.4%)
0.02198
1.298
0.99973
9
2,395,769
3:24
0.02903 (+2.2%)
0.02265
1.282
0.99972
Bracket
3
2,387,731
0:54
1.9702 ± 0.0315
1.0861 ± 0.0352
1.815
0.99833
5
2,387,725
1:40
1.9806 ± 0.0328
1.1317 ± 0.0272
1.750
0.99830
Appendix
Table 26: Spline degree on single-scale L-SplineFNO (FiLM conditioning) across the benchmark cases. Heated-duct and FDS run under region mode basis . Errors, L2/L1 and R2 are mean ± s.d. over seeds 42, 43 and 44 where a cell has several (train/test split and initialisation vary with the seed); training time is the seed mean. For the multi-variable cases (Aircraft, FDS) R2 is the mean of the per-variable scores. Degrees 7 and 9 are a single seed (42), so their absolute errors are not comparable with the seed means above them; the bracketed figure is the change against degree 3 at the same seed.
Figure 19: Gradient comparisons across autograd vector Jacobian products across the entire model, Decoder queries and the analytic spline across the six benchmark cases.
Autograd VJP
Decoder autograd
Analytic spline
(coord. sensitivity)
(pointwise ∇u , latent fixed)
Geometry
u
N
ms
MiB
ms
MiB
ms
MiB
Heat sink
T
16 384
13.3
235.8
4.7
125.8
2.5
39.9
Bracket
σ
8 192
13.7
86.3
4.1
49.8
2.2
20.5
Aircraft
P
131 072
58.9
423.1
14.2
248.4
3.9
83.4
FDS
ΔT
98 496
31.7
362.6
10.2
237.4
3.7
84.3
Appendix
Table 27: L-SplineFNO spatial-derivative cost: time and additional peak memory. Decoder autograd and the analytic spline compute the same pointwise gradient of field u with the encoded latent grid fixed. Autograd VJP measures full-model coordinate sensitivity. Timings cover only derivative evaluation and unit conversion, excluding encoding and prediction, with a fresh forward state per call. Values are medians of 20 runs after 3 warm-ups on an A100 80GB (batch 1, FP32, TF32 off). Memory is the extra peak CUDA allocation above the prepared state, including any saved autograd graph, not total GPU memory.
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
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.
Lucas Gerken Starepravo, Georgios Fourtakas, Steven Lind +3
School of Engineering, The University of Manchester, Manchester, UK · School of Engineering, Cardiff University, Cardiff, UK
Neural operators (NOs) are designed to learn maps between infinite-dimensional function spaces. We propose a novel reframing of their use. By introducing an auxiliary base-space, any finite-dimensional function can be viewed as an operator acting by composition on functions of the base-space. Through a range of benchmarks on analytic functions of increasing complexity and dimensionality, we demonstrate that NOs can match or outperform standard multilayer perceptrons and Kolmogorov--Arnold Networks in accuracy while requiring significantly fewer parameters and training time. As a real-world application, we apply a two-dimensional Tensorized Fourier Neural Operator (TFNO) to the nuclear chart, learning a correction to state-of-the-art nuclear mass models as a partially observed residual field. A TFNO ensemble reaches a held-out root-mean-square error of 198.2 keV, placing it among the best recent neural-network approaches while retaining high parameter efficiency and short training times. More broadly, these results introduce NOs as a scalable framework for finite-dimensional function interpolation, from analytic benchmarks to structured scientific data.
Department of Physics, CCTP and ITCP, University of Crete, 71303, Greece · Theoretical Physics Department, CERN, Geneva, Switzerland · Physik-Institut, Universität Zürich, 8057 Zürich, Switzerland