Learnable Hodge-Laplacian Driven Mixture of Experts for Steady-PDE Surrogates on Coarse Unstructured Meshes

Stavros Nousias; Akis Nousias; Konstantinos Gkrispanis

cs.LGOct 14, 2026

Learnable Hodge-Laplacian Driven Mixture of Experts for Steady-PDE Surrogates on Coarse Unstructured Meshes

Authors: Stavros Nousias, Akis Nousias, Konstantinos Gkrispanis

Organizations: Technical University of Munich · K3Y Labs

Abstract

Current approaches for steady partial differential equations surrogate models on unstructured meshes apply a single mixing mechanism to the whole solution. This paper proposes an architecture that partitions the solution among vanilla transformer experts, each focused on a different spectral scale. The partition is built from a discrete Hodge Laplacian $L_0=\star_0^{-1} d_0^{\top}\star_1 d_0$ , that serves as a fitting layer measuring high-frequency content. A gating mechanism takes the mesh as input and produces the per-vertex weights that recombine the expert outputs. AirFRANS is a steady Reynolds-averaged Navier--Stokes (RANS) airfoil benchmark. We evaluate our approach in a reduced remeshed version of the benchmark with each unstructured mesh counting 9.6K faces and 4.9K vertices. We compare our approach with the Transolver on the same dataset and dataset slit. Our approach reduces relative $L_2$ error by 31% and interior mean squared error by 47%.

Figures & tables

model	params	mesh use	$\text{rel-}L_{2}\downarrow$	$R^{2}\uparrow$	Volume $\downarrow$	Surf $\downarrow$	$C_{L}\downarrow$	$\rho_{L}\uparrow$
Transolver [ 17 ]	$6.00$ M	none	0.0436	0.9953	0.00265	0.00143	0.0103	0.99989
PointNet [ 43 ]	$0.08$ M	none	0.0689	0.9904	0.00451	0.00440	0.0149	0.99965
Pointwise MLP	$0.02$ M	none	0.1743	0.9404	0.03576	0.03969	0.0860	0.99338
PSD-star HodgeFormer [ 21 ]	$0.47$ M	operator	0.0800	0.9885	0.01028	0.01048	0.0681	0.99883
Hodge-routed MoE (ours)	$3.02$ M	router	0.0300	0.9983	0.00140	0.00164	0.0225	0.99988

Table 1: Comparison on the official AirFRANS

800/200

split. Every row is evaluated on the same

\approx 4.9

k-vertex working meshes and evaluation metrics. Baselines use their full input-feature arms; training and checkpoint-selection differences are specified in section 4.1 . Mesh use names how mesh structure reaches the prediction: none for a point-set model, operator when a learned Hodge Laplacian computes the fields, router when it computes only the mixing weights.

C_{L}

is the relative error of the mesh-native lift over the

200

test cases and

\rho_{L}

its Spearman correlation with the reference.

Figure 2: Pressure $p$ (top) and turbulent viscosity $\nu_{t}$ (bottom) on the reference colour scales; left to right, the CFD reference, (a) our model, (b) Transolver, (c) PointNet, and (d) the best PSD-star HodgeFormer. Values beneath each model panel are mean relative $L_{2}$ errors over the $200$ test cases. The models differ more on $\nu_{t}$ than on $p$ ; the reference panel shows fine-scale structure in the $\nu_{t}$ wake.

Figure 3: One test case at the reported checkpoint (epoch $380$ ). Columns are the reference, the prediction and the absolute error; rows are $u_{x}$ , $u_{y}$ , $p$ and $\nu_{t}$ . All panels are on the model’s own working mesh in physical units, with the reference and the prediction sharing each field’s colour scale. The per-field relative $L_{2}$ for this case is $0.0047$ , $0.0260$ , $0.0281$ and $0.0676$ . Test-set means are listed in table S7 ; the panel illustrates one case rather than a test-set average.

symbol	type	meaning	symbol	type	meaning
The discretisation
$\mathcal{M}$	$(V,E,F)$	the working mesh as a simplicial complex	$e=(i,j)$	index	an edge with endpoints $i$ , $j$
$V,E,F$	sets	vertices, edges, faces	$i,j,k,m,s$	indices	elements $i,j$ ; expert $k$ ; atom $m$ ; stream $s$
$N,N_{e},N_{f}$	integers	$\|V\|,\|E\|,\|F\|$ ; $N\approx 4900$	$c$	index	output channel in $y_{vc}$
$v$	index	a vertex
Input features

Table S1: Notation. Principal symbols grouped by their role. Array shapes omit the batch axis; local uses in the HodgeFormer preliminaries and transformer equations are identified explicitly.

optimiser	AdamW [ 41 ] , weight decay $10^{-4}$
schedule	OneCycleLR per batch [ 42 ] , with $\max_{\text{lr}}=10^{-3}$ , pct_start $=0.3$ , div $=25$ , final_div $=10^{4}$ , that is $4\times 10^{-5}\to 10^{-3}$ over the first $30\%$ and cosine decay thereafter
epochs / batch	$400$ / one whole mesh
precision	bfloat16 autocast; incidence products and the whole analyzer forced to fp32
dropout	$0.1$ throughout
initialisation	Xavier uniform on every parameter of dimension greater than one
checkpoints	every $20$ epochs; per-epoch metrics appended to metrics.jsonl

Table S2: Base training configuration. The VEF variant adds a

200

-epoch continuation at learning rate

2\times 10^{-4}

part	symbol	value
inputs	$(c_{v},c_{e},c_{f})$	$(9,21,13)$ channels on vertices, edges and faces
embedding	$d$	$192$ ; $W_{1}:9\to 192$ , GELU, $W_{2}:192\to 192$ ( eq. 4 )
experts	$K$	$8$
	$d_{x}$	$88$ , with $W^{k}_{\text{in}}:(d+2n_{F})=208\to 88$ and $\operatorname{proj}_{k}:88\to 192$
	$L_{x}$ , $h_{x}$	$2$ pre-norm blocks, $8$ heads of dimension $11$ , feed-forward $4d_{x}=352$
	$n_{F}$ , $\sigma$	$8$ Fourier features; $\sigma=(1.0,1.7,2.8,4.8,8.0,13.5,22.6,38.0)$

Table S3: Architecture. Shared settings of the reported vertex-only model and the VEF extension. Dropout is

0.1

everywhere and every parameter of dimension greater than one is Xavier-uniform initialised.

Figure S1: Training the reported model. (a) Training objective, test-set field loss, and auxiliary band loss, with the learning rate on the right axis. The training objective includes auxiliary terms absent from the test-set field loss. (b) The

20

-checkpoint test sweep, with the selected checkpoint marked. Selection uses the unrounded mean relative

L_{2}

error.

name used in this paper	what it is	params
mesh-free predictors , where mesh use is none
Plain transformer	$4$ -layer pre-norm softmax transformer, width $64$ , global attention over all vertices; the Hodge operator deleted	$0.31$ M
windowed	attention restricted to the geodesic-ring support used by the operator	$0.31$ M
$+$ locality bias	global attention $+$ a learned soft distance bias on the logits	$0.31$ M
Wide plain transformer	the same at width $192$ (FFN $768$ , $12$ heads); the width ablation	$2.72$ M
$+$ lift loss	$+$ a differentiable $C_{L}$ contour integral in the objective	$2.72$ M

Table S4: The model family. Every trained variant in this paper, grouped by how mesh structure reaches the prediction. Router means the operator produces only gate logits; bands means the mesh enters only through the supervisory decomposition of the target; operator means the mesh is in the prediction path; none means a point-set model.

arm	channels	params	$\text{rel-}L_{2}\downarrow$	$R^{2}\uparrow$	Volume $\downarrow$	Surf $\downarrow$	$C_{L}\downarrow$	$\rho_{L}\uparrow$
pi	pos $+$ inlet	$5{,}995{,}332$	0.0439	0.9954	0.00270	0.00140	0.0078	0.99988
pia	$+$ area	$5{,}995{,}844$	0.0592	0.9922	0.00418	0.00282	0.0115	0.99980
pisn	$+$ SDF $+$ normals	$5{,}996{,}868$	0.0440	0.9952	0.00286	0.00137	0.0098	0.99977
pisna	$+$ SDF $+$ normals $+$ area	$5{,}997{,}380$	0.0436	0.9953	0.00265	0.00143	0.0103	0.99989

Table S5: Transolver by input arm , own-grid metrics on the official

800/200

split.

model	mesh use	params	$\text{rel-}L_{2}\downarrow$	$R^{2}\uparrow$	Volume $\downarrow$	Surf $\downarrow$	$C_{L}\downarrow$
Hodge-routed MoE	router	$3.02$ M	0.0300	0.9983	0.00140	0.00164	0.0225
wide experts	router	$5.68$ M	0.0341	0.9979	0.00163	0.00219	0.0286
VEF router	router	$3.11$ M	0.0359	0.9978	0.00182	0.00180	0.0277
Learned-band MoE	bands	$2.75$ M	0.0345	0.9979	0.00165	0.00192	0.0265
Fixed-band MoE	bands	$2.75$ M	0.0351	0.9978	0.00178	0.00184	0.0280
$+$ wall-surface branch	bands	$3.00$ M	0.0378	0.9973	0.00196	0.00227	0.0247

Table S6: Complete comparison. Official AirFRANS

800/200

split, own-grid evaluation on the

{\approx}4.9

k-vertex mesh, the same test meshes and evaluation metrics for every row. Objectives differ across loss ablations; training and selection protocols are specified in section 4.1 . Models are described in section S5 . Best per column in bold .

model	$u_{x}$	$u_{y}$	$p$	$\nu_{t}$
Hodge-routed MoE	0.0056	0.0259	0.0255	0.0629
wide experts	0.0072	0.0319	0.0337	0.0638
VEF router	0.0067	0.0348	0.0380	0.0642
Learned-band MoE	0.0068	0.0288	0.0365	0.0658
Fixed-band MoE	0.0068	0.0312	0.0346	0.0679
$+$ wall-surface branch	0.0064	0.0293	0.0342	0.0812

Table S7: Per-field relative

L_{2}

for selected models. Between the learned-band MoE and the Hodge-routed MoE, the three smooth fields

u_{x}

u_{y}

and

p

decrease by

17\%

10\%

and

30\%

, and

\nu_{t}

decreases by

4\%

	Volume (interior)				Surf (wall)
model	$u_{x}$	$u_{y}$	$p$	$\nu_{t}$	$u_{x}$	$u_{y}$	$p$	$\nu_{t}$
Hodge-routed MoE	0.00011	0.00052	0.00031	0.00466	0.00073	0.00174	0.00248	0.00160
VEF router	0.00019	0.00096	0.00077	0.00534	0.00079	0.00197	0.00272	0.00173
Learned-band MoE	0.00019	0.00063	0.00072	0.00506	0.00087	0.00195	0.00307	0.00181
Fixed-band MoE	0.00019	0.00086	0.00069	0.00539	0.00090	0.00198	0.00296	0.00152
Wide plain transformer	0.00017	0.00109	0.00074	0.00859	0.00062	0.00191	0.00313	0.00201

Table S8: Per-field Volume and Surf (AirFRANS-native normalised MSE, best checkpoint). For the Hodge-routed MoE,

\nu_{t}

accounts for

83\%

of the summed Volume error over the four fields.

model	checkpoints	best $\text{rel-}L_{2}$ (epoch)	best Volume (epoch)
Hodge-routed MoE	20	0.0300 (380)	0.00139 (400)
wide experts	20	0.0341 (400)	0.00163 (400)
VEF router †	30	0.0359 (560)	0.00182 (560)
Learned-band MoE	20	0.0345 (380)	0.00165 (380)
Fixed-band MoE	20	0.0351 (400)	0.00178 (380)
Plain transformer	20	0.0544 (320)	0.00339 (400)

Table S9: Checkpoint sweep. Both headline metrics with the epoch at which each is lowest. † the VEF router was trained for

400

epochs plus a

200

-epoch continuation at

\mathrm{lr}=2\times 10^{-4}

and swept at

30

checkpoints rather than

20

	Hodge-routed MoE ( $V$ )	VEF router ( $VEF$ )	change
operators	$L_{0}$	$L_{0}$ , $L_{1}^{\text{down}}\!+\!L_{1}^{\text{up}}$ , $L_{2}^{\text{down}}$
supports	vertex/edge $64/64$	vertex/edge/face $64/32/1$
params	$3{,}023{,}962$	$3{,}106{,}490$	$+82{,}528$ ( $+2.7\%$ )
epochs / checkpoints	$400$ / $20$	$400+200$ / $30$
best checkpoint	$380$	$560$
$\text{rel-}L_{2}\downarrow$	0.0300	0.0359	$+19.7\%$

Table S10:

V

versus

VEF

routing. The reported controller updates vertices alone (

V

); VEF updates vertices, edges, and faces. VEF has

2.7\%

more parameters and was trained

50\%

longer, with a different edge support cap. Errors are evaluated at the checkpoint selected by mean relative

L_{2}

over

200

test cases. Percentage changes use unrounded metrics.

Figure S2: The nine telescoping atoms of the analyzer for the illustrated pressure field. Atom $0$ is the coarsest remainder; atoms $1$ – $7$ are differences between smoothing levels, and atom $8$ is the finest residual. The coarsest and finest components account for approximately $81\%$ and $14\%$ of the summed squared atom norms. Their sum reconstructs the field algebraically, subject to floating-point rounding ( eq. 14 ).

Figure S3: The analyzer and assignment for the illustrated case. Edge conductance $w_{e}$ ; the $8\times 9$ assignment $P$ ( eq. 15 ); the eight assigned targets on a shared symmetric colour range; and normalized squared component norms before and after assignment. Global Xavier initialization overwrites the constructor’s preferred assignment logits.

Figure S4: The gate $g(v)$ for the illustrated case. Top: dominant expert at each vertex, per-vertex gate entropy, and mean gate mass per expert against the uniform $1/K$ line (vertical axis restricted to $0.120$ – $0.130$ ). Bottom: the eight weight fields $g_{k}(v)$ . These panels describe routing for the displayed mesh rather than averages over the test set.

Figure S5: The fixed reference decomposition of the pressure field with the ladder $(2,8,32)$ on the unweighted adjacency, for comparison with fig. S2 . The normalized squared component norms for this illustrated field are $[0.65,0.06,0.06,0.23]$ ; the components overlap.

PDE Surrogate Modeling

PDE surrogate modeling builds fast approximations of partial differential equation solutions for simulation, prediction, and design.

See all 325 papers

Structure-Preserving Neural Surrogates with Tractable Uncertainty Quantification
Jun 10, 2026Handi Zhang, Adrienne M. Propp, Brooks Kinch +2Uncertainty Quantification PDE Surrogate Modeling
Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, USA · Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA · Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA, USA +1
Accelerating PDE Surrogates via RL-Guided Mesh Optimization
Mar 2, 2026Yang Meng, Ruoxi Jiang, Zhuokai Zhao +3PDE Surrogate Modeling Data-Efficient Learning
UChicago · Fudan University · UAlbany
When Attention Beats Fourier: Multi-Scale Transformers for PDE Solving on Irregular Domains
May 8, 2026Brandon Yee, Pairie Koh, Jack Rodriguez +1Neural PDE Solvers PDE Solving
Physics Lab, Yee Collins Research Group

Jun 10, 2026J/K move · Enter open · S save

Structure-Preserving Neural Surrogates with Tractable Uncertainty Quantification

Handi Zhang, Adrienne M. Propp, Brooks Kinch, Houman Owhadi+1

Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, USA · Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA · Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA, USA+1

Uncertainty Quantification PDE Surrogate Modeling Reduced-Order Modeling Finite Element Method Gaussian Process Regression

Open paper AskPDF arXiv

Recent advances in scientific machine learning provide a means of near-real-time solution to partial differential equations (PDEs), but lack the theoretical underpinnings of conventional simulators that support contemporary verification and validation. In this work, we construct data-driven reduced-order models that serve as structure-preserving, real-time surrogates. Remarkably, the exterior calculus that imposes physical conservation structure also exposes topological structure that we use to build a Gaussian process (GP) representation of uncertainty in state-flux relationships, ultimately yielding a Dirichlet-to-Neumann map for quantities of interest with closed-form expressions for posterior uncertainty. We specifically propose structure-preserving

H(\mathrm{div})

L^2

subspaces of conventional Raviart--Thomas and

dgP_0

elements prescribed by a lightweight transformer. Reduced-order dynamics consistent with this subspace are learned by posing a conservation law in which a GP describes the fluxes between volumes. This work hinges on a novel interface between mixed FEM spaces and GP regression; when training is posed as the optimal recovery problem (ORP), the resulting GP regression can be written as an optimization problem with equality constraints that impose a conservation structure, amenable to a fast Schur-complement training strategy. The trained model can then be solved in real time with closed-form estimators for boundary fluxes driven by prescribed Dirichlet data. The paper includes RKHS posterior error bounds for linear functionals to support uncertainty quantification, as well as numerical experiments demonstrating the accuracy of the posterior distribution as a surrogate for error estimation.

Abstract

PDE Surrogate Modeling

Explore similar work

Structure-Preserving Neural Surrogates with Tractable Uncertainty Quantification

Accelerating PDE Surrogates via RL-Guided Mesh Optimization

When Attention Beats Fourier: Multi-Scale Transformers for PDE Solving on Irregular Domains