Topology-Preserving Neural Operator Learning via Hodge Decomposition
Dongzhe Zheng 1 Tao Zhong 1 Christine Allen-Blanchette 1
arXiv:2605.13834v1 [cs.LG] 13 May 2026
Abstract
et al., 2003; Arnold et al., 2006; 2010); a linear-algebraic primer for readers is provided in Appendix A. Within this framework, the exterior derivative d : Ωk → Ωk+1 , the codifferential δ : Ωk → Ωk−1 , and the Hodge star operator ∗ : Ωk → Ωn−k uniformly characterize gradient, divergence, curl, and Laplacian operators. Many PDEs L can be written as A(u; g, κ, ∂M, f ) = 0, where u ∈ k Ωk (M) is a multi-order differential form field, κ is a material property tensor, f is a source term, g is the metric tensor, ∂M is the manifold boundary, and A is obtained by combining d, δ, and ∗ with κ. From the perspective of operator learning, finding a numerical solution to A can be viewed as learning a continuous operator G : (f, u|∂M , κ) 7→ u that is reusable across meshes and geometries.
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.
In purely Euclidean domains, neural operator methods have achieved significant successes: Fourier Neural Operators realize global convolution through low-rank spectral kernels (Li et al., 2020a), DeepONet approximates operators via dual-branch encodings (Lu et al., 2021), and PINNs embed PDE residuals directly in losses (Raissi et al., 2019). These approaches leverage regular grids and fast spectral transforms for resolution-independent operator approximation (Kovachki et al., 2023). However, many critical applications involve physical fields on Riemannian manifolds with boundaries, curvature, and non-trivial topology—including aerodynamic fields on vehicle surfaces, geophysical fields on spherical manifolds, and biological fields on organ geometries. Such quantities correspond to differential forms whose evolution is jointly constrained by cohomological structure and Riemannian metric, making them sensitive to discretization choices. Constructing neural operators on general Riemannian manifolds that are both resolutionindependent and structure-preserving constitutes the core problem this work addresses.
1. Introduction Problem Background A wide range of continuum physics models (e.g., fluid mechanics, elasticity, electromagnetic fields, and reaction-diffusion systems) can be uniformly represented as partial differential operator equations on Riemannian manifolds with a boundary (Kovachki et al., 2023). Given a finite-dimensional Riemannian manifold (M, g) with a boundary, physical fields on M are represented as differential forms of various orders: 0-forms correspond to scalar fields such as temperature or potential energy, 1-forms correspond to flux-type covector fields such as mass flow rate or current density, and 2-forms correspond to fluxes through surface elements or vorticity (Hirani, 2003; Desbrun
Research Problem and Challenges On a Riemannian manifold (M, g) with a boundary and non-trivial topology, the temporal evolution of physical fields is simultaneously constrained by two fundamentally different types of structural constraints: topological and geometric. The inherent tension between preserving global structure and resolving local dynamics constitutes the core design trade-off in the design of our method.
1
Princeton University. Correspondence to: Dongzhe Zheng <[email protected]>, Christine Allen-Blanchette <[email protected]>. Proceedings of the 43 rd International Conference on Machine Learning, Seoul, South Korea. PMLR 306, 2026. Copyright 2026 by the author(s).
Topological constraints stem from Hodge theory: the ker1
Topology-Preserving Neural Operator Learning via Hodge Decomposition
nel of the Hodge Laplacian ∆k = dδ + δd gives harmonic forms isomorphic to the k-th cohomology group, encoding global invariants such as circulation and net flux that must be explicitly preserved (Bhatia et al., 2012; Lim, 2020). Geometric and material constraints from the Riemannian metric g and material tensor κ govern high-frequency dynamics, diffusion anisotropy, and fine-scale structures such as boundary layers. The Hodge decomposition uniquely separates each k-form into gradient-type, curl-type, and harmonic components, orthogonally decoupling local differential structure from global conservation (Bhatia et al., 2012; Lim, 2020). Discrete exterior calculus and finite element exterior calculus preserve this structure exactly on simplicial complexes (Hirani, 2003; Desbrun et al., 2003; Arnold et al., 2006; 2010): vertices, edges and faces carry discrete 0-, 1-, 2-forms, and discrete operators maintain d2 = 0, δ 2 = 0, and Hodge decomposition.
(2) A spectral–geometric dual-branch design that separates topology-constrained global structure from geometry-driven local dynamics, enabling complementary and stable approximation of physical operators; (3) Empirical results showing improved accuracy over existing neural operator methods on complex manifold geometries, with exact preservation of cohomological invariants. Our result reveals that Hodge orthogonality gives operator learning on manifolds an additive approximation property, enabling geometry-driven dynamics to complement topological structure while correctly completing spectral energy.
2. Related Work Local Methods Based on Graph and Geometric Deep Learning Graph-based approaches treat meshes as graphs, employing message passing or gauge equivariant convolution to approximate PDE-induced local coupling (Bronstein et al., 2017; Cohen et al., 2019; Weiler et al., 2021). However, local aggregation mechanisms exhibit structural bottlenecks in modeling long-range dependencies: oversmoothing and over-squashing hinder networks from capturing global topological structure determined by the Hodge Laplacian kernel (Bhatia et al., 2012; Li et al., 2018; Xu et al., 2018; Oono & Suzuki, 2019; Cai & Wang, 2020; Lim, 2020; Alon & Yahav, 2020; Wang et al., 2025). Moreover, standard GNNs lack explicit encoding of differential complexes and higher-order forms, leaving algebraic identities as soft constraints rather than architectural inductive biases.
Extending neural operators (Li et al., 2020a; Lu et al., 2021; Kovachki et al., 2023) to manifold settings reveals fundamental tensions. Intrinsic geometric methods based on geodesic or tangent bundle convolution (Bronstein et al., 2017) preserve manifold structure. However, they require geometry-adaptive kernels, incurring prohibitive overhead on large meshes and struggling with high-frequency patterns. Extrinsic spectral methods leverage FFT on Euclidean grids for efficient global convolution (Li et al., 2020a; Serrano et al., 2023), yet remain agnostic to cohomological and boundary topology, with topological invariants only softly penalized rather than architecturally enforced. Graph-based methods rely on message passing or attention (Bronstein et al., 2017), suffering from over-smoothing or quadratic complexity, while neglecting higher-order simplicial adjacencies essential for cohomological structure (Li et al., 2018; Alon & Yahav, 2020; Wang et al., 2025). These limitations indicate that embedding the differential complex (d, δ, ∆k ) as architectural inductive biases while efficiently capturing high-frequency dynamics governed by metric g and material tensor κ remains open. This raises a natural question: how can operator learning on discrete meshes jointly address higher-order differential form structure across varying geometries while avoiding the efficiency–expressiveness trade-offs and topological blind spots of current approaches?
Neural Operators and Spectral Methods on Manifolds Neural operators such as FNO and DeepONet have achieved significant progress on Euclidean domains by learning function space mappings (Raissi et al., 2019; Li et al., 2020a; Lu et al., 2021; Kovachki et al., 2023). When extending to manifolds, extrinsic embedding or background grid methods struggle to preserve intrinsic metrics and flux conservation at the discrete level (Serrano et al., 2023). Recent works attempt intrinsic operators via Laplace–Beltrami eigenbases or implicit neural fields (Serrano et al., 2023; Chen et al., 2024; Liu et al., 2025), but these primarily target scalar fields without systematically incorporating de Rham complex structure, leaving harmonic components and topological invariants implicitly entangled.
Overview of This Work This paper proposes the Hodge Spectral Duality (HSD) framework for neural operator learning on oriented simplicial complexes, transforming PDE solving into structured learning on higher-order graphs with a dual-branch architecture coupled through Lie–Trotter type operator splitting (Hairer et al., 2006; Blanes et al., 2024). Our main contributions are: (1) A structure-aware neural operator framework on simplicial complexes that incorporates discrete exterior calculus as an algebraic inductive bias, ensuring physically consistent operator learning;
Higher-Order Graphs, Discrete Exterior Calculus, and Topological Deep Learning Discrete exterior calculus (DEC) and finite element exterior calculus provide rigorous frameworks for preserving algebraic and cohomological structure on simplicial complexes (Hirani, 2003; Desbrun et al., 2003; Arnold et al., 2006; 2010; Bhatia et al., 2012; Lim, 2020), including weighted Laplacian variants (Yadokoro & Bhattacharya, 2023). Topological 2
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 1. Overview of the Hodge Spectral Duality (HSD) architecture. The HSD architecture separates operator learning into a spectral Base branch (bottom) for global topology and an ambient Fiber branch (top) for high-frequency geometry. A commutator module and orthogonal projection integrate these components, ensuring strict preservation of topological invariants on manifolds.
deep learning leverages this theory for higher-order feature learning through simplicial neural networks (Papillon et al., 2023; Zia et al., 2024; Papamarkou et al., 2024; Isufi et al., 2025; Ebli et al., 2020; Chen et al., 2022; Hajij et al., 2022). However, existing works mostly focus on classification or finite-step interpolation tasks, lacking continuous operator mapping capabilities, and few explicitly separate topologydominated from metric-dominated components.
differential forms. Denote ω k ∈ C k (K, R),
f k ∈ C k (K, R)
as the discrete unknown field and right-hand side. Given discretization Ak : C k (K, R) → C k (K, R) of continuous operator Ak , the steady-state equation is Ak ω ⋆k = f k ,
ω ⋆k = G k (f k ),
(1)
where G k : C k (K, R) → C k (K, R) is the true solution operator (for time-dependent problems, the evolution operator over a given interval). Our neural operator Gθk approximates G k directly on C k (K, R).
3. Method: Hodge Spectral Duality Operator This section constructs the Hodge Spectral Duality neural operator on simplicial complexes, treating signals as discrete physical fields: 0-forms on nodes (scalar potentials like temperature, pressure), 1-forms on edges (flows like velocity, current), and 2-forms on faces (fluxes like magnetic flux, vorticity).
Discrete exterior calculus gives the k-th order Hodge–de Rham Laplacian Lk : C k (K, R) → C k (K, R), assembled from boundary operators, discrete exterior derivative dk , codifferential δk , and Hodge star ∗k (Appendix B). Construction can use topological ML libraries TopoX/TopoNetX (Hajij et al., 2022; 2024). Offline, we solve the sparse eigenvalue problem
The approach uses Hodge Decomposition under orthogonal projection via Lie-Trotter operator splitting to decouple fields into global topological and local geometric modes, with targeted neural components for each: (1) Global Topology (Base Space): Low-frequency harmonic forms captured efficiently in the spectral domain, avoiding costly spatial long-range computations (Section 3.2); (2) Local Geometry (Ambient Fiber Space): High-frequency gradient and curl components processed via spatial convolution, exploiting local dependencies to avoid full-graph redundancy (Section 3.3). Continuous physical models and PDE operators appear in the problem description; discrete exterior calculus, Hodge Laplacian, and tangent bundle definitions are in Appendix B; complexity analysis and implementation details are in Appendix F.
Lk Ψk = Ψk Λk ,
(2)
truncating to mk eigenpairs (all harmonic modes plus lowest-frequency non-harmonic modes), yielding orthogonalized spectral basis Φk ∈ RNk ×mk . This defines base k space Vbase = span(Φk ) and orthogonal complement fiber k space Vfiber under Hodge inner product ⟨·, ·⟩∗k . Projection operators and decomposition details are in Appendix C. Fields ω k are decomposed via Hodge orthogonal projection: base space components (low-dimensional spectral coefficients) and fiber components (high-frequency, metricdominated local structures). Two complementary branches model these:
3.1. Discrete Operator Learning and Hodge Spectral Decomposition
k k Gθk = Gbase,θ + Gfiber,θ .
Let (M, g) be a compact oriented Riemannian manifold with boundary, K an oriented simplicial complex approximating it, and C k (K, R) the space of k-th order discrete
(3)
k k Here Gbase,θ : C k (K, R) → Vbase learns topologydominated low-frequency response in truncated Hodge spectral domain, preserving cohomological information and
3
Topology-Preserving Neural Operator Learning via Hodge Decomposition (ℓ+1)
k k conservation laws; Gfiber,θ : C k (K, R) → Vfiber captures metric-related high-frequency corrections via tangent bundle embedding (Appendices B, C).
3.2. Base Space Branch: Spectral Domain Coefficient Learning
3.3. Fiber Branch: Metric-Dominated Correction on Tangent Bundle
k The base space branch Gbase,θ operates within the truncated k spectral subspace Vbase : project discrete fields to Hodge spectral domain, perform physically constrained nonlinear mapping in this low-dimensional space, reconstruct to base space, achieving resolution-independent operator approximation while preserving topological structure.
k The Fiber branch Gfiber,θ captures local high-frequency dynamics dominated by the Riemannian metric g and material tensor κ (anisotropic diffusion, boundary layers) without disrupting global topology encoded by the base branch. It models residuals in an auxiliary Euclidean spectral domain and constrains outputs to the base space complement via orthogonal projection. Mathematical consistency and Reach condition constraints are detailed in Appendix D.
(ℓ)
At layer ℓ, the current field ω k ∈ C k (K, R) yields spectral coefficients via Hodge inner product: (ℓ)
(ℓ)
mk ck = Φ⊤ k ∗k ω k ∈ R
(ℓ)
The layer output reconstructs via ω k,base = Φk c̃k (complete derivation in equation (33), Appendix C), achieving k good Vbase approximation under Hodge inner product and providing a topologically consistent low-frequency anchor for the fiber branch.
(4)
We introduce structure-preserving operators between discrete form space C k (K, RCℓ ) and an auxiliary Euclidean grid Ωaux . The lift operator ι : C k (K) → L2 (Ωaux ) extends discrete cochains to ambient tensor fields via Whitney forms and kernel density estimation. The pullback operator R : L2 (Ωaux ) → C k (K) maps ambient fields back through trilinear interpolation and Whitney projection, forming an adjoint pair with ι under discrete Hodge inner product.
where Φk is the truncated spectral basis from equation (2) (ℓ) (formalized in Appendix C). The coefficient vector ck maintains dimension mk across mesh resolutions, encoding harmonic and low-frequency non-harmonic modes. To embed discrete differential structure, the DEC operators dk and δk are pre-projected onto the truncated basis, yield(k) (k) ing spectral derivative matrices Md and Mδ (matrix forms in Appendix B). The branch constructs combined features: (ℓ) (ℓ) (k) (ℓ) (k) (ℓ) qk = concat ck , Md ck , Mδ ck .
Spectral convolution in ambient space uses the FNO architecture. At layer ℓ, FFT on the auxiliary grid captures metric-related high-frequency correlations: (ℓ) (ℓ) (ℓ) ω̃ k,geom = R ◦ F −1 Rloc F ◦ ι ω k , (6)
(1) (ℓ)
For a 1-form input (e.g., velocity ω 1 ), Md c1 generates the spectral representation under exterior derivative (2-form (ℓ) vorticity/flux). Thus qk provides explicit (k + 1)-order (curl-type) and (k − 1)-order (divergence-type) derivative information while updating only k-th order coefficients.
(ℓ)
where Rloc is the learnable frequency-domain spectral kernel and F is FFT on Ωaux . This handles local geometric details via global convolution on a fixed Cartesian grid, avoiding costly anisotropic manifold convolutions. To ensure geometric corrections preserve global conservation, we introduce orthogonal projection under the discrete Hodge inner product ⟨α, β⟩H k = α⊤ Hk β. Let Πkbase k project onto Vbase ; the Fiber output is constrained to orthogk onal complement Vfiber : (ℓ) (ℓ+1) ω k,fiber = I − Πkbase ω̃ k,geom . (7)
To capture quadratic nonlinear coupling (e.g., convection terms u · ∇u), we design a gated operator gMLPk with content and gating branches: (ℓ) (ℓ) (ℓ) (ℓ) c̃k = Wout ϕ(Wg qk ) ⊙ (Wc qk ) + ck , (5) where Wg , Wc , Wout are learnable projections, ϕ is SiLU activation, and ⊙ is Hadamard product. This gating introduces multiplicative inductive bias in spectral space, approx(k) (k) imating nonlinear mode mixing from Md and Mδ .
This constraint ensures the Fiber branch corrects only highfrequency metric-dominated degrees of freedom; any lowfrequency artifacts or conservation-violating modes are eliminated by projection, guaranteeing global topological invariance during cross-scale evolution.
To preserve harmonic invariants in ker Lk , hard constraints are imposed on zero-eigenvalue modes after spectral update. k Let IH be the harmonic mode indices; a diagonal projection (ℓ) k PH replaces corresponding components of c̃k with origi(ℓ) nal ck , strictly preserving cohomology classes and global flux invariance per layer. Definition of PkH and its relation to Betti numbers bk are in Appendix C.
3.4. Commutator Error and Spectral-Geometric Coupling The commutator [AkTopo , AkGeom ] = AkTopo AkGeom − AkGeom AkTopo ̸= 0 implies operator non-commutativity: the 4
Topology-Preserving Neural Operator Learning via Hodge Decomposition Table 1. Comparison of main experimental results. All experiments are conducted under comparable parameter counts (∼207k–310k). Best results are marked in bold, and second-best results are underlined. Enstrophy fidelity is not computed for the scalar field task.
MSE↓
Grad Fid↑
Enst Fid↑
Energy Fid↑
Spec Fid↑
Sβ0 ↑
IoU↑
External Aerodynamics
GNO MGN DeepONet Geo-FNO FNO-3D HSD (Ours)
231k 247k 239k 253k 227k 207k
5.40 × 10−2 5.63 × 10−2 2.54 × 10−2 3.20 × 10−2 1.80 × 10−2 1.08 × 10−2
0.8203 0.8599 0.9448 0.9263 0.9663 0.9742
0.6025 0.4671 0.3698 0.2131 0.6638 0.7658
0.4328 0.5138 0.7431 0.6438 0.7779 0.8906
0.3004 0.4209 0.6731 0.5728 0.7110 0.8423
0.4941 0.3982 0.3625 0.2261 0.5584 0.6112
0.2910 0.2720 0.0908 0.0551 0.3010 0.3398
GNO MGN DeepONet Geo-FNO FNO-3D HSD (Ours)
231k 247k 239k 253k 227k 212k
4.33 × 10−2 3.31 × 10−3 2.89 × 10−4 8.97 × 10−4 8.51 × 10−4 1.84 × 10−4
0.5337 0.9713 0.9978 0.9930 0.9940 0.9982
0.2392 0.7016 0.8246 0.8666 0.8731 0.9444
0.3708 0.8692 0.9646 0.9044 0.9054 0.9662
0.2373 0.7561 0.9468 0.8555 0.8660 0.9492
0.2491 0.5400 0.7877 0.6703 0.7041 0.8176
0.1415 0.5653 0.7834 0.7085 0.7236 0.8110
GNO MGN DeepONet Geo-FNO FNO-3D HSD (Ours)
286k 303k 273k 284k 309k 246k
3.14 × 10−2 5.23 × 10−4 2.48 × 10−3 1.55 × 10−3 5.55 × 10−4 3.56 × 10−4
0.5056 0.7485 0.5700 0.7182 0.8006 0.8578
— — — — — —
0.3082 0.5232 0.2074 0.4104 0.6365 0.6968
0.0788 0.7712 0.7044 0.7278 0.9079 0.9115
0.3795 0.5457 0.3871 0.4631 0.6721 0.7829
0.3011 0.5758 0.4433 0.5494 0.7515 0.8131
order of applying topological and geometric operators yields different results, causing systematic splitting residuals. We (ℓ) introduce correction operator Cθ constrained to Vkfiber via k (I − Πbase ), ensuring corrections act only on high-frequency components (Appendix E).
tion on unstructured meshes. MeshGraphNets (MGN) (Pfaff et al., 2020) uses an encoder-processor-decoder architecture with MLP encoding and iterative message passing, designed for physical simulations. DeepONet (Lu et al., 2021) adopts branch-trunk decomposition: the branch encodes input function values at sensor locations, the trunk encodes query coordinates, with outputs computed via inner product. Fourier Neural Operator (FNO) (Li et al., 2020a) parameterizes kernels in the frequency domain for global dependencies; for irregular meshes, spectral convolution is performed after trilinear scattering to a Cartesian grid. Geometry-Adaptive FNO (Geo-FNO) (Li et al., 2020a) learns diffeomorphic mappings from physical coordinates to a uniform latent domain to handle geometric irregularities. All baselines are reproduced from official implementations or torch geometric/neuraloperator libraries; hyperparameters are in Appendix G.
Interaction features couple geometric lift with spectral derivatives: (ℓ) (ℓ) (k) (ℓ) (k) (ℓ) z(ℓ) = ι ω k ⊕ concat ck , Md ck , Mδ ck , (8) where (ck , Md ck , Mδ ck ) recover first-order derivatives (ℓ) (dk ω, δk ω) in discrete sense. Implementing Cθ as a lightweight MLP: (ℓ+1)
ωk
Topological Fidelity
Params
Magnetostatics
Physics Consistency
Model
Toroidal Transport
Standard Accuracy
Task
(ℓ)
(ℓ)
= Gbase ω k
h i (ℓ) (ℓ) (ℓ) + (I − Πkbase ) Gfiber ω k + Cθ z(ℓ) .
(9) Near-zero initialization allows gradual learning of commutator-dominated coupling from a decoupled state.
4.2. Evaluation Metrics We construct a multi-dimensional evaluation framework encompassing accuracy, physical conservation, and topological consistency. Mean Squared Error (MSE) measures pointP (i) (i) wise distance: MSE = N −1 i ∥upred − ugt ∥2 . Gradient Fidelity (Grad Fid) evaluates gradient consistency between predicted and ground truth fields. Spectral Fidelity (Spec Fid) computes weighted relative error in the Hodge Laplacian spectral domain: Fidspec = exp(−α∥Λ−1/2 (ĉpred − ĉgt )∥/∥Λ−1/2 ĉgt ∥), where ĉ denotes Hodge spectral coefficients and Λ is the eigenvalue matrix.
4. Experiments We evaluate HSD on three tasks spanning geometric complexity, topological connectivity, and dynamic evolution: flow field reconstruction, magnetostatic field solving in multiply-connected domains, and transport processes with periodic topology; results are summarized in Table 1. 4.1. Baseline Models We compare HSD against five mainstream neural operator methods. Graph Neural Operator (GNO) (Li et al., 2020b) defines kernel integration on graphs via radius neighborhood message aggregation to approximate continuous convolu-
Physical conservation metrics use discrete exterior differential operators. REnstrophy Fidelity (Enst Fid) compares enstrophy Eω = |ω|2 (where ω = ∇ × u is vorticity) to 5
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 2. Visualization of velocity vector field predictions for the External Aerodynamics task. Columns correspond to different models, with the top row showing predictions and the bottom row showing corresponding pointwise absolute errors.
detect non-physical vorticity dissipation; applicable only to vector field tasks. Energy R Fidelity (Energy Fid) compares Dirichlet energy ED = |∇u|2 , reflecting velocity field energy distribution for vector tasks and concentration gradient intensity for scalar tasks.
continuous operators on surfaces (Bronstein et al., 2017), with sparse linear solvers handling anisotropic diffusion problems (Desbrun et al., 2003). Table 1 (upper) shows results for this task. MSE of 1.08 × 10−2 (40% reduction vs FNO-3D). On enstrophy fidelity (0.7658) and spectral fidelity (0.8423)—key metrics for turbulent microstructures and energy distributions, HSD effectively avoids over-smoothing and captures highfrequency vortex features. Geo-FNO’s coordinate deformation mechanism struggles with complex industrial geometries under constrained parameter scales, losing subtle physical features. Additional visualizations are in Appendix H.2.
Topological consistency is evaluated via β0 Score and level set IoU. The β0 Score measures connected component consistency across thresholds: Sβ0 = Eλ [exp(−|β0 (uλpred ) − β0 (uλgt )|/ max(β0 (uλgt ), 1))], quantifying capture of independent physical features. IoU measures the isosurface geometric overlap for the spatial accuracy of high-response regions. For vector fields, topological metrics are computed on vorticity for stable features.
4.4. Magnetostatic Fields in Multiply-Connected Domains (Magnetostatics)
4.3. Incompressible Flow Reconstruction on Complex Geometries (External Aerodynamics)
This task simulates magnetostatic field problems excited by source distributions in three-dimensional space, with the physics described by the steady-state form of Maxwell’s equations (Jackson, 2021). The magnetic flux density B is modeled as a linear superposition of scalar potential gradients and non-local harmonic fields induced by topological defects B = −∇ϕ + Bharm , where the scalar potential satisfies the Poisson equation ∆ϕ = ρm (Bhatia et al., 2012). The operator learning task is defined as mapping from the source scalar field (magnetic charge density ρm ) to the divergence-free magnetic flux vector field B (Kovachki et al., 2023).
This task simulates incompressible viscous fluid flow defined on curved manifolds (Marsden & Ratiu, 2013), with the objective of reconstructing the velocity field from vorticity distributions. The mathematical core involves solving the Laplace-Beltrami equation ∆M ψ = ζ on the manifold, followed by velocity field reconstruction through orthogonal gradient decomposition u = ∇⊥ ψ + uharm (Bhatia et al., 2012), strictly satisfying the divergence-free constraint ∇ · u = 0. The operator learning task is defined as mapping from the scalar vorticity field ω on the manifold surface to the velocity vector field u in the tangent space (Kovachki et al., 2023).
The computational domain is a three-dimensional bounding box containing spherical shell obstacles, forming a multiplyconnected region. This geometric structure is discretized using unstructured tetrahedral meshes (Hang, 2015), sampled to 3000 nodes. The presence of spherical shell obstacles makes the computational domain non-simply connected, and the magnetic field must contain global flux components induced by geometric cavities in addition to the irrotational component generated by local sources (Nakahara, 2003). Numerical solutions employ finite element methods to discretize the governing equations, combined with multipole expansion techniques to construct global harmonic bases
Geometric data are sourced from high-fidelity automotive triangular meshes in the DrivAerNet++ dataset (Elrefaie et al., 2024), sampled to 3000 nodes. This scenario exhibits high geometric complexity, containing non-smooth features, sharp edges, and regions with curvature variations. On closed surfaces, fluid dynamics are constrained by global geometric properties (Arnold et al., 2010), where the velocity field solution depends not only on local vorticity distributions but also on global circulation constraints determined by the surface genus (Hatcher, 2002). Numerical solutions employ graph discrete differential operators to approximate 6
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 3. Slice visualization of magnetic vector field at z = 0 plane for the Magnetostatics task. Columns correspond to different models, with the top row showing predictions and the bottom row showing corresponding errors, with the color scale: black→red→yellow→white indicates Figure 4. Comparison of Magnetostatincreasing error. ics field predictions.
tiotemporal modeling capabilities beyond local receptive fields (Gu & Yau, 2008). Numerical solutions employ a hybrid strategy, with the advection term using finite volume methods with flux limiters in upwind schemes (Versteeg & Malalasekera, 2007), and the diffusion term using implicit time integration schemes (Hairer et al., 1993).
Figure 5. Topological contours: level set connectivity (threshold at 50% of range). The β0 values for each model are: GT=3, HSD=3, GNO=2,DeepONet=0.
Table 1 (lower) shows results for this task. MGN achieves MSE of 5.23 × 10−4 , demonstrating message passing effectiveness for local temporal evolution. HSD achieves MSE of 3.56 × 10−4 (36% reduction vs FNO-3D), with Energy Fidelity (0.6968 vs 0.6365) and β0 Score (0.7829 vs 0.6721) significantly outperforming FNO-3D, indicating stable preservation of energy dissipation and topological structure under long-time integration while avoiding nonphysical artifacts (Li et al., 2020a). Additional visualizations are in Appendix H.3.
satisfying boundary conditions (Jin, 2015). Table 1 (middle) shows results for this task. Since magnetostatics is dominated by global topological obstacles rather than high-frequency dynamics, DeepONet achieves lower MSE and optimal spectral fidelity via global fitting. HSD maintains comprehensive advantages: MSE of 1.84 × 10−4 (36% reduction vs DeepONet) and enstrophy fidelity of 0.9444, precisely identifying local flux concentration structures that DeepONet statistically smooths out, preserving energy distribution and topological features (Karniadakis et al., 2021). Additional visualizations are in Appendix H.1.
4.6. Spectral Bias Analysis Figure 7 shows the predicted field energy distribution on Hodge Laplacian eigenmodes. DeepONet, Geo-FNO, and GNO exhibit rapid energy decay with increasing eigenfrequency, falling below ground truth, indicating highfrequency filtering. HSD’s spectrum closely matches the ground truth in high-frequency regions, overcoming spectral bias. The Fiber branch effectively compensates the Base branch’s high-frequency deficiencies, validating the spectral-geometric duality design.
4.5. Advection-Diffusion Dynamics (Toroidal Transport) This task simulates time-varying transport processes of scalar fields driven by a given velocity field (LeVeque, 2002), governed by the unsteady advection-diffusion equation ∂u/∂t + ∇ · (vu) = ν∆u. The system couples hyperbolic advective transport mechanisms with parabolic diffusive dissipation mechanisms, exhibiting significant spatiotemporal multiscale characteristics (Quarteroni, 2009). The operator learning task is defined as mapping from the initial scalar concentration distribution u0 to the temporal evolution trajectory of the scalar field u(t).
4.7. Training Efficiency and Computational Complexity Despite introducing spectral decomposition preprocessing, HSD maintains significant training efficiency advantages. The key is decoupling offline geometric encoding from online learning: Hodge Laplacian eigendecomposition executes once per mesh, taking approximately 57 seconds on the most complex tetrahedral mesh (Magnetostatics, ∼20k elements) and only seconds on surface meshes.
The computational domain is a torus embedded in threedimensional space, topologically possessing non-zero genus (Genus=1) with two independent non-contractible closed loop paths (Hatcher, 2002), sampled to 3000 nodes. The multiply-connected nature of the torus gives the scalar field transport process a re-entrant characteristic, where physical quantities must maintain continuity after traversing around the torus, requiring the model to possess long-range spa-
Online complexity is O(N k) for spectral projection and O(N log N ) for FFT, versus message-passing MGN’s O(N |E|). Including all preprocessing, HSD’s total training 7
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 6. Comparison of initial conditions, final ground truth field, and scalar field predictions from each model. Columns correspond to different models, with the top row showing predicted fields and the bottom row showing corresponding errors relative to Ground Truth, with the same color scale as Figure 3.
(a) Magnetostatics
(b) External Aerodynamics
(c) Toroidal Transport
Figure 7. Spectral energy decay analysis of predicted fields across the three tasks. The horizontal axis represents eigenfrequency λ, and the vertical axis represents spectral coefficient energy |ck |2 . Table 3. Impact of spectral truncation number k on model performance (MSE). Percentages in parentheses indicate change relative to k = 64.
Table 2. Ablation study results for core components (MSE). Percentages in parentheses indicate performance degradation relative to complete HSD. Model Variant HSD (Full) w/o Cθ w/o Πbase FNO-3D
Magnetostatics
Ext. Aero.
Toroidal Trans.
Modes k
Magnetostatics
Ext. Aero.
Toroidal Trans.
1.84 × 10−4 2.18 × 10−4 (+18%) 2.20 × 10−4 (+20%) 8.51 × 10−4 (+363%)
1.08 × 10−2 1.17 × 10−2 (+8%) 1.45 × 10−2 (+34%) 1.80 × 10−2 (+67%)
3.56 × 10−4 3.79 × 10−4 (+6%) 3.72 × 10−4 (+4%) 5.55 × 10−4 (+56%)
k = 64 k = 128 k = 256
1.84 × 10−4 1.64 × 10−4 (↓11%) 1.58 × 10−4 (↓14%)
1.08 × 10−2 9.32 × 10−3 (↓14%) 8.99 × 10−3 (↓17%)
3.56 × 10−4 2.80 × 10−4 (↓21%) 2.73 × 10−4 (↓23%)
ization under remeshing: when the inference mesh density increases from 3000 to 7000 nodes, its error varies by at most 30%, whereas all baselines suffer at least a 10× larger error amplification, indicating HSD learns the underlying physical operator rather than mesh-specific mappings.
time is 56× faster than MGN in External Aerodynamics (33s vs 1865s) and only 5% of MGN in Magnetostatics (215s vs 3983s), while remaining comparable to purely Euclidean FNO-3D. Complete comparisons are in Appendix Table 6.
5. Discussion, Limitations and Scope 4.8. Ablation Studies
This method relies on sparse spectral decomposition of the discrete Hodge Laplacian, requiring fixed geometry, isomorphic/isometric deformations, or minor perturbations to amortize eigendecomposition costs to offline precomputation. The framework is thus currently suited for Eulerianperspective simulations on manifolds of three dimensions or less. For scenarios requiring per-step mesh topology reconstruction, future work will incorporate iso-spectral deformation or Functional Maps theory to enable low-cost spectral basis transfer across time-varying geometry and non-isometric deformations without repeated eigendecomposition. Additionally, the low-pass characteristics of ambient mollification accommodate high-gradient continuous structures (e.g., boundary layers) but not strong discontinuities such as shock waves, restricting applicability to shockfree regimes. Despite these constraints, the proposed archi-
Table 2 analyzes core component contributions. Removing orthogonal projection causes spectral convolution to introduce non-physical low-frequency noise, with the largest degradation on geometrically complex domains. Removing the commutator MLP most impacts multiply-connected domains, confirming that topological-geometric operator non-commutativity requires compensation via Cθ . Table 3 shows that increasing spectral modes improves performance with diminishing returns, validating the spectralgeometric duality: the Base branch needs only a few lowfrequency modes for global topology, while the Fiber branch efficiently captures high-frequency details without costly large eigenbases. On External Aerodynamics, HSD exhibits stable general8
Topology-Preserving Neural Operator Learning via Hodge Decomposition
tecture successfully addresses the difficulty of traditional neural operators in simultaneously capturing high-frequency geometric details and global conservation laws on complex manifolds. This provides empirical evidence for an intrinsic Hodge-orthogonal additive structure underlying operator learning on manifolds.
hedral cnn. In International conference on Machine learning, pp. 1321–1330. PMLR, 2019. Desbrun, M., Hirani, A. N., and Marsden, J. E. Discrete exterior calculus for variational problems in computer vision and graphics. In 42nd IEEE international conference on decision and control (IEEE cat. no. 03CH37475), volume 5, pp. 4902–4907. IEEE, 2003.
Impact Statement
Ebli, S., Defferrard, M., and Spreemann, G. Simplicial neural networks. arXiv preprint arXiv:2010.03633, 2020.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
Elrefaie, M., Morar, F., Dai, A., and Ahmed, F. Drivaernet++: A large-scale multimodal car dataset with computational fluid dynamics simulations and deep learning benchmarks. Advances in Neural Information Processing Systems, 37:499–536, 2024.
References Alon, U. and Yahav, E. On the bottleneck of graph neural networks and its practical implications. arXiv preprint arXiv:2006.05205, 2020.
Federer, H. Curvature measures. Transactions of the American Mathematical Society, 93(3):418–491, 1959. Gu, X. and Yau, S. Computational Conformal Geometry. Advanced lectures in mathematics. International Press, 2008. ISBN 9781571461711.
Arnold, D., Falk, R., and Winther, R. Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American mathematical society, 47(2): 281–354, 2010.
Hairer, E., Nørsett, S., and Wanner, G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems. Springer, 1993. ISBN 9783540604525.
Arnold, D. N., Falk, R. S., and Winther, R. Finite element exterior calculus, homological techniques, and applications. Acta numerica, 15:1–155, 2006.
Hairer, E., Lubich, C., and Wanner, G. Structure-preserving algorithms for ordinary differential equations. Geometric numerical integration, 31, 2006.
Bhatia, H., Norgard, G., Pascucci, V., and Bremer, P.-T. The helmholtz-hodge decomposition—a survey. IEEE Transactions on visualization and computer graphics, 19 (8):1386–1404, 2012.
Hajij, M., Zamzmi, G., Papamarkou, T., Miolane, N., Guzmán-Sáenz, A., Ramamurthy, K. N., Birdal, T., Dey, T. K., Mukherjee, S., Samaga, S. N., et al. Topological deep learning: Going beyond graph data. arXiv preprint arXiv:2206.00606, 2022.
Blanes, S., Casas, F., and Murua, A. Splitting methods for differential equations. arXiv preprint arXiv:2401.01722, 2024. Bronstein, M. M., Bruna, J., LeCun, Y., Szlam, A., and Vandergheynst, P. Geometric deep learning: going beyond euclidean data. IEEE Signal Processing Magazine, 34(4): 18–42, 2017.
Hajij, M., Papillon, M., Frantzen, F., Agerberg, J., AlJabea, I., Ballester, R., Battiloro, C., Bernárdez, G., Birdal, T., Brent, A., et al. Topox: a suite of python packages for machine learning on topological domains. Journal of Machine Learning Research, 25(374):1–8, 2024.
Cai, C. and Wang, Y. A note on over-smoothing for graph neural networks. arXiv preprint arXiv:2006.13318, 2020.
Hang, S. Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Softw, 41(2):11, 2015.
Chen, G., Liu, X., Meng, Q., Chen, L., Liu, C., and Li, Y. Learning neural operators on riemannian manifolds. National Science Open, 3(6):20240001, 2024.
Hatcher, A. Algebraic Topology. Cambridge University Press, 2002.
Chen, Y., Gel, Y. R., and Poor, H. V. Bscnets: Block simplicial complex neural networks. In Proceedings of the aaai conference on artificial intelligence, volume 36, pp. 6333–6341, 2022.
Hirani, A. N. Discrete exterior calculus. California Institute of Technology, 2003. Isufi, E., Leus, G., Beferull-Lozano, B., Barbarossa, S., and Di Lorenzo, P. Topological signal processing and learning: Recent advances and future challenges. Signal Processing, pp. 109930, 2025.
Cohen, T., Weiler, M., Kicanaoglu, B., and Welling, M. Gauge equivariant convolutional networks and the icosa9
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Jackson, J. D. Classical electrodynamics. John Wiley & Sons, 2021.
Oono, K. and Suzuki, T. Graph neural networks exponentially lose expressive power for node classification. arXiv preprint arXiv:1905.10947, 2019.
Jin, J. The Finite Element Method in Electromagnetics. IEEE Press. Wiley, 2015. ISBN 9781118842027.
Papamarkou, T., Birdal, T., Bronstein, M., Carlsson, G., Curry, J., Gao, Y., Hajij, M., Kwitt, R., Lio, P., Di Lorenzo, P., et al. Position: Topological deep learning is the new frontier for relational learning. Proceedings of machine learning research, 235:39529, 2024.
Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L. Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021.
Papillon, M., Sanborn, S., Hajij, M., and Miolane, N. Architectures of topological deep learning: A survey of message-passing topological neural networks. arXiv preprint arXiv:2304.10031, 2023.
Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1–97, 2023.
Pfaff, T., Fortunato, M., Sanchez-Gonzalez, A., and Battaglia, P. Learning mesh-based simulation with graph networks. In International conference on learning representations, 2020.
LeVeque, R. J. Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, 2002.
Quarteroni, A. Numerical models for differential problems. Springer, 2009.
Li, Q., Han, Z., and Wu, X.-M. Deeper insights into graph convolutional networks for semi-supervised learning. In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physicsinformed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020a.
Serrano, L., Le Boudec, L., Kassaı̈ Koupaı̈, A., Wang, T. X., Yin, Y., Vittaut, J.-N., and Gallinari, P. Operator learning with neural fields: Tackling pdes on general geometries. Advances in Neural Information Processing Systems, 36: 70581–70611, 2023.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485, 2020b.
Versteeg, H. and Malalasekera, W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education Limited, 2007. ISBN 9780131274983.
Lim, L.-H. Hodge laplacians on graphs. Siam Review, 62 (3):685–715, 2020.
Wang, K., Yang, Y., Saha, I., and Allen-Blanchette, C. Resolving oversmoothing with opinion dissensus. arXiv preprint arXiv:2501.19089, 2025.
Liu, S., Liu, H., Zhang, T., and Liu, X. Ms-iuffno: Multiscale implicit u-net enhanced factorized fourier neural operator for solving geometric pdes. Computer Methods in Applied Mechanics and Engineering, 437:117761, 2025.
Weiler, M., Forré, P., Verlinde, E., and Welling, M. Coordinate independent convolutional networks–isometry and gauge equivariant convolutions on riemannian manifolds. arXiv preprint arXiv:2106.06020, 2021.
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218–229, 2021.
Xu, K., Hu, W., Leskovec, J., and Jegelka, S. How powerful are graph neural networks? arXiv preprint arXiv:1810.00826, 2018.
Marsden, J. E. and Ratiu, T. S. Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems, volume 17. Springer Science & Business Media, 2013.
Yadokoro, S. and Bhattacharya, S. Weighted combinatorial laplacian and its application to coverage repair in sensor networks. arXiv preprint arXiv:2312.04825, 2023. Zia, A., Khamis, A., Nichols, J., Tayab, U. B., Hayder, Z., Rolland, V., Stone, E., and Petersson, L. Topological deep learning: a review of an emerging paradigm. Artificial Intelligence Review, 57(4):77, 2024.
Nakahara, M. Geometry, Topology and Physics, Second Edition. Graduate student series in physics. Taylor & Francis, 2003. ISBN 9780750306065. 10
Topology-Preserving Neural Operator Learning via Hodge Decomposition
A. Notation Reference Table and A Primer on Algebraic Topology Symbol (M, g)
Type Continuous geometry
∂M K Kk , Nk = |Kk | C k (K, R)
Continuous geometry Discrete geometry Discrete geometry Space
ω k ∈ C k (K, R) f k ∈ C k (K, R) Ak
State Data Continuous operator
Ak : C k → C k Gk Gθk d, δ ∆k = dδ + δd Bk
Discrete operator Solution operator Neural operator Continuous operator Continuous operator DEC operator
dk = B ⊤ k+1 ∗k
DEC operator DEC operator
⟨α, β⟩∗k = α⊤ ∗k β δk = ∗−1 k−1 Bk ∗k Lk = dk−1 δk + δk+1 dk bk = dim ker Lk
Inner product DEC operator DEC operator Topological quantity
Lk Ψk = Ψk Λk ψk,i λk,i Ψk Λk mk
Spectral decomposition Mode Eigenvalue Matrix Matrix Truncation dimension
Φk ∈ RNk ×mk k Vbase = span(Φk )
Basis Subspace
k k Vfiber = (Vbase )⊥∗k
Subspace
Πkbase = Φk Φ⊤ k ∗k ω k,base = Πkbase ω k ω k,fiber = (I − Πkbase )ω k AkTopo AkGeom k k Gbase : C k → Vbase k k k Gfiber : C → Vfiber k k Gbase,θ , Gfiber,θ
Projection Component Component Continuous operator Continuous operator Solution operator Solution operator Neural operator
(ℓ)
ωk (ℓ) (ℓ) ck = Φ ⊤ k ∗k ω k (k) ⊤ Md = Φk+1 ∗k+1 dk Φk (k) Mδ = Φ⊤ k−1 ∗k−1 δk Φk (ℓ) qk MLPk (ℓ) c̃k k IH PkH
Network state Spectral coefficients Matrix Matrix Features Network Spectral coefficients Index set Projection
Meaning Compact oriented Riemannian manifold with boundary and its Riemannian metric. Boundary of manifold M. Oriented simplicial complex approximating (M, g). Set of k-dimensional simplices and its cardinality. Space of k-th order discrete differential forms (k-cochains), carrying discrete degrees of freedom of k-th order physical quantities. Discrete representation of unknown k-th order physical field. Discrete representation of right-hand side/source term or external force. Control operator acting on k-th order differential forms in continuous PDE. Discretization of Ak on C k (K, R). (True) solution operator satisfying ω ⋆k = G k (f k ). Neural operator with parameters θ, used to approximate G k . Exterior derivative d and codifferential δ. k-th order Hodge–de Rham Laplacian. Oriented boundary matrix between k-dimensional and (k − 1)dimensional simplices. Discrete exterior derivative C k (K, R) → C k+1 (K, R). k-th order discrete Hodge star operator (mass matrix), inducing discrete inner product. Discrete Hodge inner product on C k (K, R). Discrete codifferential C k (K, R) → C k−1 (K, R). k-th order discrete Hodge Laplacian. k-th Betti number, corresponding to the dimension of harmonic k-form space. Eigendecomposition of discrete Hodge Laplacian. i-th eigenvector of Lk (k-th order spectral mode). i-th eigenvalue of Lk . Matrix composed of all eigenvectors ψk,i . Diagonal eigenvalue matrix diag(λk,1 , . . . , λk,Nk ). Number of modes corresponding to the smallest eigenvalues retained in truncated spectrum (including all harmonic modes). Truncated spectral basis [ψk,1 , . . . , ψk,mk ], orthogonalized under ∗k . Base space spectral subspace containing harmonic modes and several low-frequency modes. k Hodge orthogonal complement of Vbase , high-frequency / local degree of freedom subspace. k Orthogonal projection operator onto Vbase . Base space (topological + low-frequency) component of ω k . Fiber (high-frequency / metric-dominated) component of ω k . Topology-dominated operator component generated by d, δ, ∆k . Geometry / material-dominated operator component depending on g, κ. Solution operator mapping input to base space subspace. Solution operator mapping input to fiber subspace. k k Learnable approximations of Gbase and Gfiber . Current field approximation at layer ℓ. (ℓ) Coordinates of ω k under truncated spectral basis. Discrete exterior derivative operator under truncated spectral basis. Discrete codifferential operator under truncated spectral basis. (ℓ) (k) (ℓ) (k) (ℓ) Spectral domain features: concat(ck , Md ck , Mδ ck ). (ℓ) Small multilayer perceptron acting on qk . Spectral coefficients after update by MLPk . Spectral index set corresponding to zero eigenvalues (harmonic modes). Diagonal projection matrix onto harmonic subspace. (Continued on next page)
11
Topology-Preserving Neural Operator Learning via Hodge Decomposition Symbol (ℓ+1) (ℓ) ω k,base = Φk c̃k ι
Type Component Lift operator
R Ωaux
Pullback operator Domain
F , F −1 (ℓ) Rloc
Operator Kernel
(ℓ)
Field Component
Cθ
(ℓ)
Correction operator
z(ℓ)
Features
λ
Scalar
τM
Geometric quantity
h, haux ϵ γ
Scale Scale Mesh quality
ω̃ k,geom (ℓ+1) (ℓ) ω k,fiber = (I − Πkbase )ω̃ k,geom
Meaning Field component reconstructed from base space branch at layer ℓ + 1. Lifts discrete forms to continuous fields on a manifold/auxiliary Euclidean domain (Whitney forms + kernel smoothing). Interpolates / projects ambient fields back to C k (K, R). Auxiliary Euclidean domain surrounding the manifold (domain of regular voxel grid). Discrete fast Fourier transform and its inverse on Ωaux . Frequency domain convolution kernel tensor for ambient FNO at layer ℓ. Geometric residual field after ambient FNO correction. k Fiber branch output after projection onto Vfiber . Learnable correction term approximating commutator error [AkTopo , AkGeom ]. (ℓ) Spectral-geometric interaction features (concatenation of ι(ω k ) and (ℓ) qk ). Global scaling/confidence coefficient for Fiber geometric residual (ResScale). Reach of manifold M, controlling geometric safety radius of ambient embedding. Simplicial complex mesh scale and auxiliary voxel grid step size. Ambient kernel smoothing bandwidth (mollification scale). Maximum condition number of element geometric matrices, characterizing mesh anisotropy.
To facilitate the understanding of the Hodge Spectral Duality framework within the machine learning community, this section reformulates the core concepts of Discrete Exterior Calculus (DEC) and Algebraic Topology using standard linear algebra and graph signal processing terminology. We emphasize exact mathematical definitions over heuristic metaphors. A.1. Data Representation: From Node Signals to Cochains In standard Graph Neural Networks (GNNs), data is typically treated as signals on vertices (X ∈ R|V |×c ). In our framework, physical fields are strictly typed by their integration domains, formalized as Discrete Differential Forms (or Cochains). • Linear Algebra Perspective: The space of k-forms, denoted as C k (K, R), is simply a vector space RNk , where Nk is the number of k-dimensional simplices (vertices, edges, faces). – 0-forms (u ∈ RN0 ): Scalars defined on vertices (e.g., Temperature). – 1-forms (v ∈ RN1 ): Scalars defined on oriented edges (e.g., Flow rate along a pipe). Note that changing edge orientation negates the value. – 2-forms (w ∈ RN2 ): Scalars defined on oriented faces (e.g., Flux through a surface element). • Insight for ML: Unlike GNNs, which learn arbitrary feature vectors, this framework enforces a strict dimension-binding inductive bias. A velocity field must be processed as a 1-form (edge signal), not a 0-form (node feature), to preserve its transformation properties under geometric deformation. A.2. The Algebraic Structure: Boundary and Coboundary The core connectivity of the mesh is encoded in the Boundary Operator ∂k . • Definition: ∂k maps a k-simplex to a linear combination of its (k − 1)-faces. • Matrix Representation: In the discrete setting, this is exactly the incidence matrix Bk ∈ RNk−1 ×Nk . For example, for an edge eij = [vi , vj ], the boundary is vj − vi . Thus, B1 is the standard edge-to-vertex incidence matrix containing only {0, 1, −1}. 12
Topology-Preserving Neural Operator Learning via Hodge Decomposition
The Exterior Derivative dk (used for Gradient, Curl) is the adjoint (transpose) of the boundary operator: ⊤ dk = Bk+1
(10)
• d0 (Gradient): Maps vertex signals to edge signals (computes differences). • d1 (Curl): Maps edge signals to face signals (sums circulation around a face). The Fundamental Property: The defining algebraic structure of a complex is ∂k−1 ◦ ∂k = 0. In matrix terms: Bk Bk+1 = 0
=⇒
dk+1 dk = 0
(11)
This identity (d2 = 0) strictly guarantees that “the Curl of a Gradient is zero” and “the Divergence of a Curl is zero” at machine precision purely via sparse matrix multiplication, without needing to learn these physics constraints. A.3. The Generalized Laplacian and Hodge Decomposition Standard GCNs utilize the Graph Laplacian L0 = D − A ≈ B1⊤ B1 . Hodge Theory generalizes this to higher dimensions: Lk = dk−1 δk + δk+1 dk | {z } | {z } Grad-Div term
(12)
Curl-Curl term
• Spectral Interpretation: Lk is a symmetric positive semi-definite matrix. Its eigenvectors provide a Fourier basis for signals on edges (k = 1) or faces (k = 2). • Hodge Decomposition: Just as any vector can be projected onto orthogonal axes, any discrete field ω ∈ RNk decomposes orthogonally into three subspaces determined by the operators above: ω = im(dk−1 ) ⊕ im(δk+1 ) ⊕ ker(Lk )
(13)
This separates the signal into Irrotational (gradient-flow), Solenoidal (divergence-free), and Harmonic components. A.4. Betti Numbers: Topological Invariants as Null Spaces Betti numbers (bk ) are often cited abstractly, but in our computational framework, they have a precise linear algebraic definition related to the Cohomology Groups. • Definition: The k-th Betti number is the dimension of the kernel (null space) of the k-th Hodge Laplacian. bk = dim(ker(Lk )) = number of zero eigenvalues of Lk
(14)
• Physical Meaning in R3 : – b0 : Number of connected components. A harmonic 0-form is constant on each component. – b1 : Number of independent non-contractible loops (e.g., flow circulating around a handle or hole). A harmonic 1-form represents a circulation that cannot be explained by a local gradient potential. – b2 : Number of enclosed voids (cavities). A harmonic 2-form represents flux trapped on a closed surface. Relevance to Operator Learning: In standard neural networks, global topological features (like b1 circulation) often get smoothed out by local message passing. By explicitly projecting onto the kernel of Lk (the harmonic subspace), our method preserves these global invariants as hard constraints, ensuring the network respects the fundamental topology of the physical domain.
B. Mathematical Foundations of Discrete Exterior Calculus and Tangent Bundle This appendix provides the strict mathematical definitions of discrete exterior calculus, Hodge spectral structure, and tangent bundle, supporting the construction of the Hodge Spectral Duality neural operator in the main text. For related theory, see (Hirani, 2003; Desbrun et al., 2003; Arnold et al., 2006; 2010). 13
Topology-Preserving Neural Operator Learning via Hodge Decomposition
B.1. Simplicial Complex and Discrete Differential Forms Let (M, g) be a compact oriented n-dimensional Riemannian manifold with boundary. To approximate (M, g) on a computer, take an oriented simplicial complex embedded in Euclidean space K = (V, E, F, . . . ), where V, E, F are the sets of vertices, edges, and faces respectively, corresponding to 0, 1, 2-dimensional simplices. In general, denote Kk as the set of k-dimensional simplices and Nk = |Kk | as its cardinality. Cases where k > n are not considered. The space of discrete k-th order differential forms C k (K, R) is defined as the real vector space of all mappings from each k-dimensional simplex σ ∈ Kk to a real number: C k (K, R) ≃ RNk , where each component corresponds to an integral quantity on a k-dimensional simplex. C 0 (K, R) corresponds to discrete scalar fields on vertices, C 1 (K, R) corresponds to line integral fluxes on oriented edges, C 2 (K, R) corresponds to area fluxes or area densities on oriented faces, and higher-order spaces C k (K, R) follow analogously. To represent adjacency and orientation relationships between simplices, for k ≥ 1, define the oriented boundary matrix Bk ∈ RNk−1 ×Nk , whose (i, j)-th component is ( 0, [Bk ]ij = ±1,
i σk−1 ̸⊂ ∂σkj , i σk−1 ⊂ ∂σkj ,
where the sign is determined by the relative orientation between simplices. Matrix Bk encodes the oriented (k − 1)dimensional boundary of each k-dimensional simplex. B.2. Discrete Exterior Derivative, Hodge Star, and Codifferential The discrete counterpart of the continuous exterior derivative operator d : Ωk (M) → Ωk+1 (M) is determined by the algebraic properties of the boundary operator. Define the discrete exterior derivative dk : C k (K, R) → C k+1 (K, R) at the matrix level as dk = B⊤ k+1 .
(15)
For k = 0, d0 gives the discrete gradient operator on graph vertices; for k = 1, d1 gives a discrete curl-type operator on edges. From Bk Bk+1 = 0 we obtain dk+1 dk = 0, reflecting the complex structure of the discrete exterior derivative. To introduce a metric-consistent inner product between discrete forms, define the discrete Hodge star operator ∗k : C k (K, R) → C k (K, R), which at the matrix level is a symmetric positive definite matrix, typically constructed through volume ratios between simplices and their dual cells or finite element mass matrices. Given ∗k , introduce the inner product on C k (K, R) ⟨α, β⟩∗k = α⊤ ∗k β,
α, β ∈ C k (K, R).
This inner product approximates at the mesh level the continuous L2 inner product Z ⟨α, β⟩ = α ∧ ∗β. M
The discrete codifferential operator δk : C k (K, R) → C k−1 (K, R) 14
(16)
Topology-Preserving Neural Operator Learning via Hodge Decomposition
is defined as the formal adjoint of dk−1 with respect to inner product (16), i.e., δk = ∗−1 k−1 Bk ∗k ,
(17)
satisfying δk2 = 0. For k = 1, δ1 corresponds to the discrete divergence operator; for k = 2, δ2 corresponds to higher-order divergence. Similar to the continuous case, the combination of dk and δk uniformly describes discrete versions of first-order differential operators such as gradient, curl, and divergence. Based on the above operators, the discrete Hodge–de Rham Laplacian is defined as Lk = dk−1 δk + δk+1 dk : C k (K, R) → C k (K, R),
(18)
whose matrix is symmetric positive semi-definite. For k = 0, if selecting standard volume metric with ∗0 = I, then L0 = B⊤ 1 B1 , which is the combinatorial graph Laplacian (unnormalized form) commonly used in graph learning. For general k, Lk generalizes the graph Laplacian to generalized Laplacian operators on higher-order cells such as edges and faces. B.3. Hodge–de Rham Decomposition and Hodge Spectrum The continuous Hodge–de Rham decomposition states that, under appropriate boundary conditions, each smooth k-form can be uniquely decomposed into three parts: gradient-type, curl/divergence-type, and harmonic-type. Discrete exterior calculus and finite element exterior calculus theory show that, on appropriate discrete shape function spaces, this decomposition still holds at the discrete level (Hirani, 2003; Desbrun et al., 2003; Arnold et al., 2006; 2010). Specifically, on C k (K, R) there is an orthogonal decomposition C k (K, R) = im dk−1 ⊕ im δk+1 ⊕ ker Lk , (19) where im dk−1 is the gradient-type component, describing parts driven by scalar potentials; im δk+1 is the curl-type or divergence-type component, describing parts driven by circulation or sources/sinks; ker Lk is the harmonic subspace, describing modes that are locally source-free but constrained by global topology. To characterize the multi-scale structure and topological modes of k-th order fields, consider the spectral decomposition of the discrete Hodge Laplacian Lk Ψk = Ψk Λk , (20) where the column vectors of Ψk ∈ RNk ×Nk form an orthogonal basis of C k (K, R), satisfying Ψ⊤ k ∗k Ψk = INk , and Λk is the diagonal eigenvalue matrix. Eigenvectors with zero eigenvalues span ker Lk , whose dimension equals the k-th Betti number bk , characterizing non-trivial topological structures such as non-contractible loops and cavities. Smaller non-zero eigenvalues correspond to large-scale modes with slow spatial variation, while larger eigenvalues correspond to localized, high-frequency modes. Given any ω k ∈ C k (K, R), its expansion under the Hodge spectral basis is ωk =
Nk X
⊤ ck,i = ψk,i ∗k ω k .
ck,i ψk,i ,
i=1
Here, i ≤ bk corresponds to harmonic modes, and i > bk corresponds to non-harmonic modes. Truncating to the first mk eigenvectors yields the spectral subspace k Vbase = span{ψk,1 , . . . , ψk,mk },
which simultaneously contains all harmonic modes and several low-frequency non-harmonic modes, used to preserve topological information and approximate global large-scale behavior in a finite-dimensional space. The error from spectral truncation is mainly concentrated in high-frequency parts, which can be corrected by local operators in the fiber branch. 15
Topology-Preserving Neural Operator Learning via Hodge Decomposition
B.4. Riemannian Geometry and Tangent Bundle Physical operators at the continuous level typically depend on local metrics, curvature, and material property tensors, which are naturally defined on the tangent bundle. The tangent bundle structure provides the foundation for constructing local operators consistent with geometry in local Euclidean coordinates. For any point p ∈ M on the manifold, the tangent space Tp M is a vector space isomorphic to Rn , whose elements are tangent vectors at p. The Riemannian metric g gives an inner product on Tp M ⟨v, w⟩g(p) = gp (v, w),
v, w ∈ Tp M,
and the corresponding volume element dvolg . The collection of tangent spaces at all points forms the tangent bundle G TM = Tp M, p∈M
which is a differentiable manifold of dimension 2n. Through local coordinate charts (Uα , φα ), φα : Uα → Rn introduces coordinates x = (x1 , . . . , xn ) within Uα , with the corresponding coordinate basis ∂ ∂ , . . . , . ∂x1 ∂xn In this basis, the metric tensor is represented as a symmetric positive definite matrix field gij (x), i.e., gp =
n X
gij (x) dxi ⊗ dxj .
i,j=1
Anisotropic diffusion tensors, stress tensors, and other quantities related to metric and material properties can be given in matrix form in local coordinates. Many first-order and second-order partial differential operators have standard expressions in local coordinates of the tangent space. Taking the scalar field u as an example, the anisotropic diffusion operator in local coordinates x can be written as n X p 1 ∂ ∂u ij ∇ · D(x)∇u(x) = p |g(x)| D (x) j , ∂x |g(x)| i,j=1 ∂xi where D(x) is a symmetric positive definite matrix field related to metric and material properties, and |g(x)| is the determinant of the metric matrix. The advection operator can be written as v(x) · ∇u(x) =
n X
v i (x)
i=1
∂u , ∂xi
where v(x) is a tangent vector field. For vector fields or higher-order form fields, the corresponding operators can be represented through combinations of covariant derivatives, exterior derivatives, and codifferentials, with coefficients also depending on local metrics and material properties. The topological properties of exterior derivative d and codifferential δ are independent of the metric, while the above diffusion and advection operators are highly sensitive to g and material property tensors. In the Hodge Spectral Duality framework, the base space branch encodes the topology-dominated parts determined by d and δ through the discrete Hodge Laplacian and its spectral structure, while the fiber branch models metric-dominated local effects through local coordinate representations on the tangent bundle and frequency domain operators, with consistency maintained through differentiable de Rham maps and orthogonal projections.
C. Spectral Operator Theory and Subspace Derivation This appendix formalizes the construction of truncated spectral subspaces, orthogonal projections, spectral derivative matrices, and harmonic projection matrices, based on the discrete exterior calculus and Hodge–de Rham spectral structure given in Appendix B, to support the simplified exposition in Sections 3.1 and 3.2 of the main text. 16
Topology-Preserving Neural Operator Learning via Hodge Decomposition
C.1. Truncated Hodge Spectral Basis and Subspace Decomposition From Appendix B, the k-th order discrete Hodge Laplacian Lk is a self-adjoint operator under the Hodge inner product ⟨α, β⟩∗k = α⊤ ∗k β,
α, β ∈ C k (K, R)
Its spectral decomposition satisfies equation (2), and the eigenvector matrix Ψk can be chosen as a Hodge orthonormal basis, i.e., Ψ⊤ k ∗k Ψk = INk . Denote the corresponding eigenvalues as 0 = λk,1 = · · · = λk,bk < λk,bk +1 ≤ · · · ≤ λk,Nk , where bk = dim ker Lk is the k-th Betti number. To obtain a low-dimensional spectral representation for operator learning, truncate to the first mk eigenvectors Lk ψk,i = λk,i ψk,i , Φk = ψk,1 , . . . , ψk,mk ∈ RNk ×mk ,
(21)
and perform a one-time orthogonalization under the Hodge inner product, so that Φ⊤ k ∗k Φk = Imk .
(22)
k Vbase = span(Φk ) ⊂ C k (K, R),
(23)
This yields the truncated spectral subspace
which necessarily contains all harmonic modes as well as several lowest-frequency non-harmonic modes. The orthogonal complement under the Hodge inner product is defined as ⊥∗k k k k Vfiber = Vbase = η ∈ C k (K, R) ⟨η, ξ⟩∗k = 0, ∀ ξ ∈ Vbase , (24) thus obtaining the approximate orthogonal decomposition k k C k (K, R) = Vbase ⊕ Vfiber .
(25)
k Based on equation (22), the Hodge orthogonal projection operator onto Vbase is
Πkbase = Φk Φ⊤ k ∗k ,
Πkfiber = I − Πkbase ,
(26)
k where Πkfiber is the projection onto the fiber subspace Vfiber . Any field ω k ∈ C k (K, R) can thus be uniquely decomposed as
ω k = ω k,base + ω k,fiber ,
ω k,base = Πkbase ω k ,
ω k,fiber = Πkfiber ω k .
(27)
Compared to the exact discrete Hodge–de Rham decomposition, equations (25)–(27) concentrate harmonic and lowfrequency modes into a finite-dimensional subspace through spectral truncation, while concentrating the remaining highfrequency degrees of freedom into the orthogonal complement. C.2. Splitting of Topological and Geometric Operators At the continuous level, the control operator Ak can be abstractly decomposed into a topology-dominated part generated by exterior derivatives and codifferentials, and a geometry-dominated part dependent on metric and material property tensors Ak = AkTopo + AkGeom ,
(28)
where AkTopo is generated by d, δ, and the Hodge–de Rham Laplacian, mainly characterizing cohomological constraints and conservation structures; AkGeom contains diffusion, advection, and source/sink terms dependent on g and material property tensors. The corresponding solution operator can be formally written as k k G k ≈ Gbase + Gfiber ,
k k Gbase : C k (K, R) → Vbase ,
k k Gfiber : C k (K, R) → Vfiber .
(29)
k k The neural operator branches Gbase,θ and Gfiber,θ in the main text equation (3) are precisely parameterized approximations of equation (29).
17
Topology-Preserving Neural Operator Learning via Hodge Decomposition
C.3. Spectral Derivative Matrices and Harmonic Projection To preserve the discrete exterior calculus structure in the truncated spectral domain, the discrete exterior derivative dk and discrete codifferential δk from DEC are projected onto the spectral basis Φk , yielding spectral derivative matrices (k)
(k)
Md = Φ⊤ k+1 ∗k+1 dk Φk ,
Mδ
= Φ⊤ k−1 ∗k−1 δk Φk .
(30)
For any field ω k , its spectral coefficients are ck = Φ⊤ k ∗k ω k . In the truncated subspace, the spectral representations of the discrete exterior derivative and codifferential approximately satisfy (k)
(k)
Φ⊤ k+1 ∗k+1 dk ω k ≈ Md ck , (k)
Φ⊤ k−1 ∗k−1 δk ω k ≈ Mδ ck ,
(k)
so Md and Mδ preserve the algebraic structure of the original differential complex in the spectral domain. The spectral features in Section 3.2 of the main text can be written as (ℓ) (ℓ) (k) (ℓ) (k) (ℓ) qk = concat ck , Md ck , Mδ ck . (31) This combined feature is then fed into gMLPk defined in the main text equation (5), using the multiplicative structure of the gating mechanism to approximate nonlinear mode coupling between physical fields. The harmonic projection is based on the zero eigenvalue structure of Lk . Denote k IH = i ∈ {1, . . . , mk } λk,i = 0 as the index set of harmonic modes in the truncated spectral basis, and define the diagonal projection matrix PkH ∈ Rmk ×mk as ( k k 1, i = j ∈ IH , PH ij = 0, otherwise, (ℓ)
then applying the harmonic hard constraint to the layer ℓ spectral update result c̃k can be written as (ℓ)
(ℓ)
(ℓ)
(ℓ)
c̃k ← c̃k + PkH ck − c̃k
,
(32) (ℓ)
which is equivalent to forcibly restoring the components at harmonic indices to the input coefficients ck , thereby exactly preserving the corresponding cohomology classes and global conserved quantities at each layer. Finally, the base space reconstruction formula is written as (ℓ+1)
(ℓ)
ω k,base = Φk c̃k ,
(33)
which together with equation (4) forms a closed loop of spectral projection and reconstruction, ensuring that the base space k branch output always lies in Vbase .
D. Consistency, Stability, and Resolution Requirements of Ambient Fiber Embedding This section analyzes the consistency of the composite operator R ◦ Aamb ◦ ιh,ϵ , consisting of discrete differential form lifting, spectral operator convolution, and pullback operations, relative to the manifold-intrinsic operator AM . Addressing the derivative discontinuity of low-order Whitney forms at element interfaces and the resolution of boundary layer features on anisotropic meshes, this section provides error bound analysis based on Reach conditions (Federer, 1959), Nyquist sampling, and anisotropic kernels, and clarifies the adaptive adjustment mechanisms for model hyperparameters. D.1. Geometric Assumptions and Discretization Setup To ensure rigor in the analysis, we first clarify the basic assumptions on geometry and discretization. Assume the manifold M has Reach τM > 0, and the simplicial complex K approximates M with Hausdorff error O(h). For mesh quality, define the condition number γK of the geometric matrix GK for each element, and let the global maximum condition number γ := supK cond(GK ) < ∞. This parameter γ characterizes the degree of mesh anisotropy, allowing elements with extreme aspect ratios, but the relevant error constants will increase with γ. 18
Topology-Preserving Neural Operator Learning via Hodge Decomposition
In the ambient embedding process, the mollification bandwidth ϵ must satisfy 0 < ϵ < cτM , where c ∈ (0, 1), to ensure uniqueness of the nearest point projection and prevent non-physical topological shortcuts. For numerical implementation, the ambient voxel grid step size haux constitutes a hard constraint on resolution. According to the Nyquist sampling theorem, haux is required to satisfy haux ≤ min(ϵ/2, c0 δres /4), where δres is the smallest characteristic scale that needs to be resolved in the physical problem. This condition indicates that the ambient grid must have sufficient resolution to support the chosen mollification scale and physical features. To mitigate aliasing problems caused by anisotropic meshes, theoretically one should adopt metric-adaptive anisotropic kernel functions, with separate settings for normal bandwidth ϵn and tangential bandwidth ϵt , satisfying ϵn ≤ min(τM /2, δres /3) and ϵt ≍ αhloc γ 1/2 . This setting ensures that high-frequency features in the normal direction are not over-smoothed while maintaining sufficient coverage in the tangential direction to suppress interpolation noise. D.2. Stability and Consistency of Whitney Extension Let Wh denote the piecewise H 1 field obtained by Whitney reconstruction from discrete forms in C k (K). Since low-order Whitney forms have derivative jumps across element interfaces, directly applying differential operators to them introduces Dirac-type singularities. The ambient mollification operator ιh,ϵ serves a regularization role here, transforming Wh into the smooth ambient field uϵ . Regarding the stability of this process, it can be shown that the gradient norm after mollification is controlled by the original discrete field, i.e., ∥∇uϵ ∥L2 (Rd ) ≤ C(γ)∥∇Wh ∥L2 (M) , where the constant C(γ) increases monotonically with mesh anisotropy. Regarding approximation consistency, for target fields u with bounded curvature and second-order regularity, the error restricted to the manifold satisfies h . (34) ∥∇(uϵ |M ) − ∇u∥L2 (M) ≤ C(γ) ϵ + ϵ Equation (34) reveals the trade-off mechanism for mollification bandwidth ϵ. The first term ϵ comes from bias introduced by mollification itself, while the second term h/ϵ reflects the residual variance from mollification failing to completely eliminate discrete mesh roughness and derivative jumps. Therefore, the theoretically optimal bandwidth should be chosen at the ϵ ≍ h1/2 level. In engineering implementation, the mollification bandwidth ϵ is not a completely independent hyperparameter, but is implicitly locked by the ambient grid resolution haux and the effective support radius of the interpolation kernel. From resolution constraints, the lower bound of the effective bandwidth is limited by voxel size, i.e., ϵeff ≥ haux . D.3. Total Error Decomposition and Adaptive Correction Combining the above geometric embedding error and spectral operator approximation error, the overall consistency of the Fiber branch can be described by total error decomposition. Let Aamb be the ambient domain FNO operator, R be the pullback operator, and assume the output is corrected by orthogonal projection (I − Πkbase ). Under the aforementioned geometric assumptions, the total error Etotal has an upper bound Etotal := ∥(I − Πkbase )(RAamb ιh,ϵ − AM )∥ ≤ C [Egeom (ϵ, h, γ) + Evox (haux ) + Espec (ξtask )] .
(35)
Here, Egeom ≍ ϵ + h/ϵ is the geometric discretization and mollification error, Evox ≍ haux is the voxelization discretization error, and Espec is the frequency domain truncation error. For cases where mesh resolution may be insufficient or aspect ratios are extremely poor, this framework achieves implicit adaptive adjustment through machine learning mechanisms. First, the learnable scalar λ in the model serves as a confidence gate. If the ambient grid resolution is insufficient to resolve high-frequency geometric features in certain regions, causing the Fiber branch output to be noisy or have large deviations, the optimization process will drive λ to decay. This allows the model to automatically degrade to a topology-preserving solution dominated by the Base branch in extreme cases, thereby ensuring numerical stability. Second, the spectral kernel weights Rloc of the ambient FNO will automatically adapt to the frequency domain truncation introduced by mollification during training. For aliased components beyond the Nyquist frequency π/haux , the network tends to learn near-zero gains, thus acting as a data-driven low-pass filter. Finally, the orthogonal projection operator (I − Πkbase ) constitutes an algebraic-level safety barrier. Even if the ambient branch introduces erroneous low-frequency modes or conservation-violating artifacts due to resolution limitations, these components will be eliminated by the projection operation. 19
Topology-Preserving Neural Operator Learning via Hodge Decomposition
The hard constraint on the harmonic subspace not only does not introduce overfitting risk, but actually constitutes a physically consistent regularization mechanism; since measurement noise typically violates conservation laws (thus falling in the im d or im δ spaces orthogonal to the harmonic kernel), enforcing topological constraints forces the model to automatically remove these non-physical residuals when minimizing the loss, thereby converging to the optimal projection of noisy data onto the conservative solution manifold. Additionally, for high-frequency discretization artifacts introduced by the voxelization process, the finite spectral bandwidth property of the ambient FNO acts as an intrinsic low-pass filter, effectively attenuating mesh-scale aliased components remaining after orthogonal projection. D.4. Efficiency and Robustness in R3 Embeddings This method primarily targets physical objects embedded in R3 , such as aerodynamic shapes and biological tissues. In these scenarios, the computational efficiency gains from Ambient FNO significantly outweigh the accuracy loss introduced by the interpolation process. The lift operator ι has a built-in controlled convolution kernel that mathematically guarantees the signal satisfies band-limited conditions before entering the ambient grid, thereby avoiding severe aliasing. Even if the ambient branch produces high-frequency noise, the projection operator I − Πbase ensures this noise is confined to the orthogonal complement space and absolutely cannot leak into the topology-dominated low-frequency subspace. This mechanism guarantees the immunity of physical conservation laws to numerical artifacts, while the remaining small high-frequency residuals are naturally smoothed through the spectral bias property of neural networks. The splatting mapping from discrete simplices to the background grid is essentially kernel density estimation. Regardless of how the manifold geometry curls or self-occludes, this process only involves local additive operations on simplices. Its computational complexity is O(Nsimplex × K 3 ), where K is the kernel width. This complexity is only linear with respect to the number of elements and completely independent of manifold curvature or topological complexity. Physically, physical fields on adjacent patches appear as a superposition or transition in ambient space, and splatting automatically handles this property without additional geometric detection costs. For extremely irregular geometries, combining sparse octree or hash grid techniques can concentrate computational resources in the manifold neighborhood, making computation proportional to manifold surface area rather than bounding box volume. D.5. Boundary Regularity and Suppression of Gibbs Artifacts Addressing the binary mask step at the manifold boundary ∂M and periodic splicing discontinuities at computational box boundaries, this section analyzes the Gibbs ringing phenomenon caused thereby and its suppression strategies. If the Fourier transform is directly applied to the binarized embedded field, spatial discontinuities will cause high-frequency coefficient decay rates to degrade from O(|k|−p ) to O(|k|−1 ), producing significant spatial domain ringing that may confuse physical high-frequency features. To address this, we introduce a processing mechanism combining geometric regularization and operator correction. First, at the geometric level, C r continuous transition regions are constructed to replace binary steps. Based on the manifold signed distance field d(x), construct a smooth soft mask mϵn (x), where ϵn is the normal bandwidth, satisfying ϵn ≲ min(τM /2, δres /3). Define the extended field uext (x) = mϵn (x)ι(u)(x) + (1 − mϵn (x))E[u](x), where E[u] is a boundary-consistent extension operator ensuring uext satisfies at least C 1 continuity at boundaries. This regularization treatment restores the Fourier coefficient decay rate to O(|k|−(r+1) ), fundamentally weakening the energy source of the Gibbs phenomenon. Meanwhile, to handle non-periodicity at computational box boundaries, smooth window functions w(x) or absorbing layers are applied at FFT domain edges to eliminate boundary artifacts introduced by periodic wrapping. Second, at the algebraic and operator level, orthogonal projection (I − Πkbase ) ensures that any artifacts falling into harmonic or low-frequency subspaces due to improper boundary handling are removed, ensuring ringing noise does not pollute global topological invariants and conservation laws. The remaining high-frequency ringing is confined to the Fiber subspace and R suppressed through introducing boundary band energy regularization terms Lring = λ |∇mϵn |2 |ωfiber |2 dx in the loss function as well as frequency domain high-frequency penalties. Under the above treatments, the total approximation error can be further decomposed into a form including boundary regularity Etot ≤ C [Egeom (ϵn , h) + Evox (haux ) + Ecut (r, ξcut ) + Ewrap (Lpml )] . (36) Here Ecut decays significantly as the transition region regularity r increases, and Ewrap decreases as the absorbing layer 20
Topology-Preserving Neural Operator Learning via Hodge Decomposition
thickness Lpml increases. This error bound indicates that by improving the regularity of the embedded field through soft masks and boundary processing, Gibbs error can be effectively controlled, validating the numerical effectiveness of this framework on bounded manifolds.
E. Theoretical Justification of the Commutator Corrector (ℓ)
This section provides the theoretical justification for the commutator correction term Cθ . We first derive the analytical form of the operator commutator [AkTopo , AkGeom ], then prove that within the Hodge Spectral Duality framework, the interaction features z(ℓ) constitute the complete derivative proxy required by this commutator. E.1. Analytical Form of the Commutator Consider k-th order differential form fields ω defined on a Riemannian manifold (M, g) with boundary. The topologydominated operator AkTopo is generated by exterior derivative complex operators (d, δ, ∆k ) and typically does not explicitly depend on local metrics; while the geometry-dominated operator AkGeom (such as anisotropic diffusion, advection) explicitly depends on position-dependent material property tensors κ(x) and Riemannian metric g(x). Using the Leibniz rule on Riemannian manifolds, the commutator of the two is typically non-zero. Taking scalar fields as an example, let ATopo = ∆ (Laplace–Beltrami operator), AGeom = Mκ (multiplication operator by κ(x)), then: [∆, Mκ ]u = ∆(κu) − κ∆u = (∆κ)u + 2 ⟨∇κ, ∇u⟩g .
(37)
The above equation shows that the residual term is driven by two parts: (1) higher-order derivatives of geometric parameters (∆κ); (2) coupling between geometric parameter gradients and physical field gradients (∇κ · ∇u). Generalizing to k-th order forms, the general form of the commutator can be written as a local differential operator Fcomm : [AkTopo , AkGeom ] ω = Fcomm ω, dω, δω, κ, dκ, g, ∇g . (38) This depends not only on zeroth-order field values ω and κ, but also on their first-order derivative information. Therefore, any network Cθ attempting to correct this spectral-geometric splitting error must explicitly or implicitly access this set of derivative information. E.2. Spectral Encoding as Derivative Proxy for Fields To capture the dependence on dω and δω in equation (38), we utilize the spectral structure of the Base branch. Recalling the (k) (k) spectral coefficients ck and spectral derivative matrices Md , Mδ defined in Section 3.2 of the main text, we have the following algebraic identities: (k) (k) dk ω k ≈ Φk+1 Md ck , δk ω k ≈ Φk−1 Mδ ck . (39) k This means that within the truncated spectral subspace Vbase , there exists a linear isomorphism between the vector group (k) (k) (ck , Md ck , Mδ ck ) and the physical field and its first-order derivatives (ω, dω, δω).
The interaction features z(ℓ) we construct in Section 3.4 explicitly concatenate this set of vectors. Therefore, for a fixed simplicial complex, z(ℓ) completely encodes the first-order differential structure of the field in the discrete sense, enabling the MLP to theoretically recover dω and δω from input features. E.3. Geometric Interaction in the Fiber Branch To capture the dependence on geometric derivatives ∇κ, ∇g in equation (38), we rely on the convolutional properties of the Fiber branch. The Fiber branch embeds discrete forms into the auxiliary Euclidean grid through the lift mapping ι and (ℓ) applies frequency domain convolution operators Rloc . Convolution operations are locally equivalent to weighted differences, (ℓ) so the output features ugeom implicitly contain the spatial variation rates (i.e., derivative information) of input parameters. In summary, the two parts of the interaction features z(ℓ) respectively provide: Base part explicitly provides algebraic derivatives of physical fields (dω, δω); Fiber part implicitly provides numerical derivatives of geometric and material property parameters (∇κ, ∇g). 21
Topology-Preserving Neural Operator Learning via Hodge Decomposition (ℓ)
Based on the Universal Approximation Theorem, the parameterized MLP Cθ can utilize this complete local information to approximate the nonlinear commutator residual Fcomm . Finally, through orthogonal projection (I − Πkbase ), we restrict this correction to the orthogonal complement space, thereby ensuring that the underlying topological conservation laws are not corrupted by approximation errors.
F. Computational Complexity Analysis F.1. Computational Complexity The Hodge Spectral Duality framework adopts an offline-online decoupling strategy, concentrating geometry-dependent spectral operations into a one-time preprocessing stage while achieving approximately linear online inference complexity with respect to mesh scale. Offline Stage: DEC Assembly and Spectral Decomposition. The offline stage constructs DEC operators on the simplicial complex K, including oriented boundary matrices Bk , discrete exterior derivatives dk , discrete Hodge stars ∗k , discrete codifferentials δk , and the discrete Hodge Laplacian Lk . These operators are highly sparse: each k-simplex is adjacent to a bounded number of (k ± 1)-simplices, yielding nnz(Lk ) = O(Nk ). Assembly involves only sparse matrix multiplications with favorable memory access patterns, achieving complexity O(Nk ). For spectral decomposition, we employ the Shift-Invert Spectral Transformation rather than direct eigenvalue computation. Instead of solving for the smallest eigenvalues of Lk directly, we solve for the largest eigenvalues of the transformed operator (Lk − σI)−1 with shift σ ≈ 0. This transformation magnifies the spectral gap between target low-frequency modes and the remainder of the spectrum, dramatically accelerating Krylov subspace convergence. Combined with sparse direct factorization of the shifted operator, the complexity for extracting mk eigenpairs becomes O(mk nnz(Lk )) ≈ O(mk Nk ), where the iteration count implicit in mk remains small due to rapid convergence. Precomputing spectral derivative matrices (k) (k) Md , Mδ and projection operators Πkbase from Φk and sparse operators dk , δk also achieves O(mk Nk ) complexity. (ℓ)
(ℓ)
Online Stage: Base Space and Fiber Branches. For the base space branch, spectral projection ck = Φ⊤ k ∗k ω k (ℓ) (ℓ) and reconstruction ω̃ k = Φk c̃k are dense matrix-vector multiplications with complexity O(Nk mk ). Since mk ≪ Nk (typically mk ∼ 64), these operations are implemented as highly parallel GEMM kernels with arithmetic intensity wellsuited for GPU acceleration. Furthermore, the exact sequence property d ◦ d = 0 and δ ◦ δ = 0 is preserved at the spectral level, allowing certain higher-order derivative compositions to be short-circuited to zero without floating-point operations. For the Fiber branch, point-voxel mapping ι interpolates signals onto a regular background grid with resolution R and total voxels V = R3 . This mapping and its inverse involve sparse interpolation with complexity O(Nk ). The 3D FNO executed on the background grid has complexity O(V log V ), where V is a fixed constant independent of mesh scale. This replaces the O(Nk ) irregular memory accesses required for per-vertex tangent space convolutions with structured FFT operations on regular grids, eliminating sparse index lookup overhead. The total complexity of a single forward pass is O(Nk mk ) + O(Nk + V log V ) ≈ O(Nk ), achieving linear scaling with mesh resolution. Compared to higher-order graph neural networks based on sparse message passing, this framework improves hardware utilization through low-dimensional spectral projection and structured FFT operations. Compared to intrinsic geometric deep learning methods, ambient space approximation decouples high-frequency geometric processing complexity from mesh resolution. F.2. Scalability and Precomputation Cost The offline-online decomposition constitutes a computational arbitrage strategy. For physical simulation tasks, geometric meshes are typically fixed or undergo isometric deformations, making the one-time preprocessing investment negligible relative to thousands of training iterations with O(Nk ) online complexity. The Shift-Invert strategy ensures that even for million-scale meshes (N ∼ 106 ), spectral decomposition completes in minutes rather than hours. Sparse direct factorization of (Lk − σI) exploits the bounded fill-in characteristic of discretizations on manifolds, and the rapid Krylov convergence induced by spectral gap magnification keeps iteration counts low regardless of mesh scale. For dynamic scenarios requiring frequent mesh topology updates, approximate spectral solvers based on multilevel algorithms 22
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Algorithm 1 Forward Pass of Hodge Spectral Duality Operator Discrete right-hand side f k , spectral basis Φk , Hodge Laplacian Lk , precomputed operators (k) (k) dk , δk , ∗k , Πkbase , Md , Mδ 2: Output: Approximate solution ω ⋆k 3: Offline Precomputation (completed once when geometry is fixed) 4: Construct boundary matrices Bk , discrete exterior derivative dk , Hodge star operator ∗k , and Hodge Laplacian Lk 5: Compute the first mk eigenpairs of Lk to obtain spectral basis Φk (k) (k) 6: Compute spectral derivative operators Md , Mδ and projection operator Πkbase 7: Online Forward Pass (0) 8: Initialize ω k (e.g., zero field or simple approximation based on f k ) 9: for ℓ = 0 to L − 1 do (ℓ) (ℓ) 10: ck ← Φ⊤ # Lift to Hodge spectral domain k ∗k ω k (ℓ) (ℓ) (k) (ℓ) (k) (ℓ) 11: qk ← concat ck , Md ck , Mδ ck # Spectral features w/ derivatives (ℓ) (ℓ) 12: c̃k ← MLPk qk (ℓ) k 13: Apply harmonic hard constraint on c̃k at IH (ℓ+1) (ℓ) 14: ω k,base ← Φk c̃k # Reconstruct base space component (ℓ) (ℓ) (ℓ+1) 15: uk ← ι ω k , ω k,base , f k , geometry # Map to ambient background grid (ℓ) (ℓ) −1 (ℓ) 16: ugeom ← F Rloc F uk # Ambient space FNO convolution (ℓ) (ℓ) 17: ω̃ k,geom ← R ugeom # Pull back to simplicial complex (ℓ) (ℓ) 18: z(ℓ) ← ι ω k ⊕ qk # Coupling features (reuse derivatives) (ℓ) (ℓ) # Nonlinear commutator correction 19: ω k,int ← Cθ z(ℓ) (ℓ) (ℓ+1) (ℓ) k 20: ω k,fiber ← I − Πbase ω̃ k,geom + ω k,int # Orthogonal projection 1: Input:
(ℓ+1)
(ℓ+1)
(ℓ+1)
21: ωk ← ω k,base + ω k,fiber 22: end for (L) 23: Return ω ⋆k ← ω k
or hierarchical matrix techniques can further reduce preprocessing overhead. Such extensions are beyond the scope of this paper, which focuses on establishing the foundational architecture.
Resolution Efficiency. A fundamental advantage of the dual-branch architecture lies in its decoupling of physical fidelity from geometric sampling density. Regardless of mesh resolution, the base space branch processes only the first mk spectral coefficients (e.g., mk = 64 or 128), confining global topological information to a fixed low-dimensional subspace whose computational cost is nearly independent of Nk . The Fiber branch handles high-frequency details through ambient space FFT with complexity O(V log V ) independent of mesh density, avoiding the global mesh refinement traditionally required to capture boundary layers or localized features. As demonstrated in our resolution robustness analysis (Section 4.8), HSD maintains consistent accuracy even when inference mesh density is significantly reduced, and models trained on coarse meshes transfer zero-shot to high-resolution meshes (7000+ vertices) with only marginal error increase. This resolution efficiency provides empirical evidence that the prohibitive computational complexity often associated with manifold PDEs is not a fundamental physical necessity. Furthermore, aggressive mesh decimation can be applied as an additional computational reduction strategy: while subsampling degrades pointwise numerical precision in high-frequency components, the dominant low-frequency spectral modes that govern global physical behavior remain well-resolved.
G. Training Details and Hyperparameter Configuration This section provides detailed information on experimental training costs, computational environment, optimization strategies, hyperparameter settings, and training procedures for each task. 23
Topology-Preserving Neural Operator Learning via Hodge Decomposition
G.1. Dataset Splitting Strategy For all three tasks, we generated 3,000 simulation samples. The dataset was split using a strict temporal partitioning strategy to prevent data leakage: the test set was first separated (20% of total data), then the remaining data was divided into training and validation sets (validation set comprising 15% of the remaining data). Detailed split statistics are shown in Table 5. Table 5. Dataset split statistics. Dataset
Proportion
Samples
Purpose
Training set
68%
2,040
Validation set
12%
360
Test set
20%
600
Gradient updates and parameter optimization Hyperparameter tuning and model selection (early stopping) Final performance evaluation and metric computation
Total
100%
3,000
—
G.2. Detailed Computational Cost Analysis Table 6 presents a detailed comparison of time and memory consumption for each model under identical hardware conditions (a single NVIDIA RTX PRO 6000 Blackwell Workstation Edition GPU). The total time for HSD includes one-time spectral basis precomputation overhead (indicated in parentheses). Table 6. Training time and computational cost comparison across tasks. Time is measured in seconds (s), and memory in megabytes (MB). HSD total time includes one-time feature basis preprocessing time (indicated in parentheses).
Magnetostatics Ext. Aerodyn. Toro. Transport Model
Params Time VRAM Time VRAM Time
VRAM
Complexity
DeepONet FNO-3D Geo-FNO HSD (Ours) GNO MGN
∼240k 6.7 ∼230k 177.0 ∼250k 49.6 ∼220k 215.5* ∼230k 477.4 ∼240k 3983.0
171.1 194.4 173.5 173.0 183.1 174.4
O(N ) O(N log N ) O(N log N ) O(N k + N log N ) O(N |E|) O(N |E|)
*
369.5 382.7 383.9 364.5 385.3 382.7
7.5 29.5 49.6 34.6* 426.5 1865.0
369.8 4.5 382.7 36.2 384.0 35.6 383.3 372.2* 385.4 195.1 383.0 1129.3
HSD time includes one-time offline spectral precomputation overhead (Magnetostatics: +57.1s, Ext. Aerodyn.: +1.6s, Toro. Transport: +3.8s).
We observe that DeepONet has the fastest training speed, but as discussed earlier, it exhibits lower physical and topological fidelity. MGN is extremely slow to train due to its explicit message passing mechanism, especially for larger mesh sizes. HSD demonstrates high efficiency in static field tasks, with slightly increased time in dynamic tasks due to the inclusion of the temporal residual correction module, but remains significantly faster than MGN and GNO. G.3. Computational Resources All experiments were implemented using the PyTorch 2.9.0 framework, specifically utilizing the NVIDIA NGC PyTorch container (Release 25.09) for optimized performance. The computations were executed on high-performance computing nodes equipped with an AMD Ryzen Threadripper 7970X 32-core CPU and dual NVIDIA RTX PRO 6000 Blackwell Workstation Edition GPUs. Each GPU provides 96 GB of VRAM, totaling 192 GB of video memory, which provides sufficient capacity for high-resolution spectral field modeling. The software environment was built upon CUDA 13.0 to ensure optimal hardware acceleration. To guarantee experimental reproducibility, we employed a Docker-based containerization strategy to maintain consistent versions of all dependency libraries and drivers. G.4. Optimization and Training Configuration All tasks employ the AdamW optimizer with a cosine annealing learning rate scheduler (CosineAnnealingLR) that gradually decays the learning rate from the initial value to zero. We designed corresponding loss functions for different physical field types: for vector field tasks (Magnetostatics and External Aerodynamics), the loss function is a weighted sum of flux loss 24
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Lflux and divergence loss Ldiv , with weights set to λflux = 1.0 and λdiv = 0.1 respectively; for the scalar field task (Toroidal Transport), we use the standard L2 relative error with an additional L1 regularization term on sparse spectral coefficients. The weight decay coefficient is set to 10−4 for vector field tasks and 10−5 for the scalar field task. G.5. Model Hyperparameters Table 7 summarizes the specific hyperparameter configurations for the HSD model and all baseline models across the three experimental tasks. To ensure fair comparison, the parameter counts for all models are controlled within a similar range (approximately 250k–300k).
Table 7. Hyperparameter configuration for each task. Hyperparameter / Model
Magnetostatics
Ext. Aero.
Toroidal Trans.
64 64 100 10−3
128 64 100 10−3
64 64 50 10−3
2 (Gated MLP) 32 43 12 4 163
2 (Gated MLP) 32 43 11 3 163
2 (Gated MLP) 32 43 11 6 163
GNO Hidden channels Projection channels Depth Neighborhood radius
84 96 5 0.2
84 96 5 0.2
120 68 3 0.15
FNO-3D Modes Hidden channels Depth Grid resolution
43 21 2 163
43 21 2 163
43 20 3 163
MGN Hidden dimension Message passing layers
58 10
58 10
72 8
[64, 64, 64] [64, 64, 64] 74
[64, 64, 64] [64, 64, 64] 74
[96, 96, 64] [68, 68, 64] 64
6 12 2
6 12 2
6 9 4
Global Settings Spectral truncation dimension k Batch size Training epochs Initial learning rate HSD (Ours) Base branch layers Base hidden dimension Fiber branch modes Fiber hidden channels Fiber depth Fiber branch grid resolution
DeepONet Branch network layers Trunk network layers Basis function dimension p Geo-FNO Modes Width Depth
Regarding hyperparameter selection, the following notes apply. In the External Aerodynamics task, the spectral truncation dimension k is increased to 128 (compared to 64 for other tasks) to accommodate the significantly higher topological complexity and surface curvature variations of DrivAerNet++ geometries. In the Toroidal Transport task, the Fiber branch depth of HSD is increased to 6 layers, as the temporal evolution characteristics of high-frequency features in the advectiondiffusion process require a deeper ambient space network to capture, whereas the static fields in Magnetostatics and External Aerodynamics tasks only require shallower network structures (4 and 3 layers, respectively). 25
Topology-Preserving Neural Operator Learning via Hodge Decomposition
G.6. Training Curves Figures 8, 9, and 10 show the training loss curves for each model across the three tasks. It can be observed that HSD exhibits stable decreasing trends in all three tasks and achieves the lowest final loss values. In contrast, FNO-3D and DeepONet show pronounced oscillations or premature plateauing in complex geometry tasks, while GNO and MGN, despite relatively stable training processes, still have higher final loss values than HSD. Notably, HSD maintains comparable training efficiency to other tasks in the External Aerodynamics task despite using a larger spectral truncation dimension.
Figure 8. Training loss curves for all models on the Magnetostatics task. The horizontal axis represents training epochs, and the vertical axis shows MSE loss on a logarithmic scale.
Figure 9. Training loss curves for all models on the External Aerodynamics task. The horizontal axis represents training epochs, and the vertical axis shows MSE loss on a logarithmic scale.
26
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 10. Training loss curves for all models on the Toroidal Transport task. The horizontal axis represents training epochs, and the vertical axis shows MSE loss on a logarithmic scale.
H. Additional Visualization Results This section provides additional visualization samples for each task to more comprehensively demonstrate the prediction performance of HSD and baseline models across different test cases. All error maps use the same color scale mapping, with black to red to yellow to white indicating increasing error magnitude. H.1. Magnetostatics Figures 11 and 12 show slice visualizations of magnetic vector fields for two additional test samples in the Magnetostatics task. It can be observed that under different geometric configurations and boundary conditions, HSD accurately captures the spatial distribution characteristics of the magnetic field, particularly maintaining low prediction errors in regions with large field strength gradients. In contrast, FNO-3D produces noticeable artifacts near complex boundaries due to its reliance on regular grid interpolation; while GNO and MGN can adapt to unstructured meshes, their field reconstruction accuracy in high-curvature regions remains inferior to HSD.
Figure 11. Slice visualization of magnetic vector field for the Magnetostatics task (Sample 1). Each column corresponds to a different model, with the top row showing predictions and the bottom row showing corresponding errors (color scale: black→red→yellow→white indicates increasing error).
27
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 12. Slice visualization of magnetic vector field for the Magnetostatics task (Sample 2). Each column corresponds to a different model, with the top row showing predictions and the bottom row showing corresponding errors. Color scale same as Figure 11.
H.2. External Aerodynamics Figures 13 and 14 show velocity vector field prediction results for two different vehicle geometries in the External Aerodynamics task. These two samples represent streamlined and bluff body designs, respectively, corresponding to distinctly different flow field characteristics. HSD demonstrates accurate prediction capability for near-wall velocity distributions in both cases, particularly in regions with geometric discontinuities such as wheel arches and side mirrors, where errors are significantly lower than those of other methods. DeepONet, due to its point-wise evaluation limitations, struggles to capture spatially correlated flow field structures; Geo-FNO, despite introducing geometric transformations, still exhibits larger errors when handling vehicle geometries with high topological complexity in the DrivAerNet++ dataset.
Figure 13. Velocity vector field prediction visualization for the External Aerodynamics task (Sample 1). Each column corresponds to a different model, with the top row showing predictions and the bottom row showing corresponding errors (color scale: black→red→yellow→white indicates increasing error).
Figure 14. Velocity vector field prediction visualization for the External Aerodynamics task (Sample 2). Each column corresponds to a different model, with the top row showing predictions and the bottom row showing corresponding errors. Color scale same as Figure 13.
Figures 15 and 16 further show surface flux distributions obtained by integrating the velocity field. As a key physical quantity for downstream engineering analysis, flux requires higher prediction accuracy. It can be seen that the flux distribution 28
Topology-Preserving Neural Operator Learning via Hodge Decomposition
predicted by HSD closely matches the Ground Truth, while baseline models exhibit pronounced spatial error accumulation effects in their flux predictions, especially in flow separation regions and wake confluence zones where errors are significantly amplified.
Figure 15. Surface flux prediction visualization for the External Aerodynamics task (Sample 1). Each column corresponds to a different model, with the top row showing predicted flux and the bottom row showing corresponding errors against Ground Truth. Color scale same as Figure 13.
Figure 16. Surface flux prediction visualization for the External Aerodynamics task (Sample 2). Each column corresponds to a different model, with the top row showing predicted flux and the bottom row showing corresponding errors. Color scale same as Figure 13.
H.3. Toroidal Transport Figures 17, 18, and 19 show three test samples with different initial conditions in the Toroidal Transport task. This task involves the advection-diffusion evolution of a scalar field on a torus, where the spatial frequency and localization of the initial field distribution directly affect the complexity of the final state. Under the low-frequency initial conditions shown in Figure 17, most models provide reasonable predictions; however, as the spatial structure complexity of the initial field increases (as in Figures 18 and 19), baseline model errors rapidly grow, manifesting as over-smoothing of high-frequency details or spurious oscillations. HSD, leveraging its spectral domain processing capability and geometry-aware architecture on fiber bundles, maintains consistently low error levels across all three samples, demonstrating its robust modeling capability for time-varying physical field evolution processes.
29
Topology-Preserving Neural Operator Learning via Hodge Decomposition
Figure 17. Scalar field prediction visualization for the Toroidal Transport task (Sample 1). From left to right: initial condition, Ground Truth, and prediction results from each model. The bottom row of each column shows corresponding errors (color scale: black→red→yellow→white indicates increasing error).
Figure 18. Scalar field prediction visualization for the Toroidal Transport task (Sample 2). From left to right: initial condition, Ground Truth, and prediction results from each model. The bottom row of each column shows corresponding errors. Color scale same as Figure 17.
Figure 19. Scalar field prediction visualization for the Toroidal Transport task (Sample 3). From left to right: initial condition, Ground Truth, and prediction results from each model. The bottom row of each column shows corresponding errors. Color scale same as Figure 17.
30