UNOP: Physics-Constrained Unsupervised Neural Operator for Long-Horizon PDE Learning on Generalized Geometries

1 Introduction

The mathematical basis for explaining physical laws is provided by partial differential equations (PDEs), which have numerous uses in energy optimization, aircraft design, and climate prediction [Raissi et al., 2019; Karniadakis et al., 2021], as well as emerging digital twin paradigms for 6G and space-ground integrated networks [Wang et al., 2023; Chen et al., 2023]. These applications involve complex spatiotemporal dynamics and generalized geometries, making traditional numerical solvers prohibitive for real-time tasks, especially in resource-constrained scenarios like mobile edge computing [Gao et al., 2019]. Inspired by the success of large-scale pre-trained models [Zeng et al., 2022], the exorbitant cost of labeled data is a major obstacle to scaling neural PDE solvers: each high-fidelity CFD/FEM simulation may need hundreds of GPU hours, rendering supervised techniques unsustainable. Consequently, shifting from label intensive supervised paradigms to physics-constrained unsupervised learning has become a necessity for scalable solvers.

Under the data-driven paradigm, three fundamental challenges emerge: (i) Label scarcity: Supervised neural operators require thousands of labeled trajectories, but generating each via CFD/FEM is expensive [Herde et al., 2024; Wu et al., 2024], creating a scalability bottleneck; (ii) Geometric constraints: Real-world domains such as vehicle surfaces, airfoils, and engine intakes [Li et al., 2023a; Li et al., 2023b; Lam et al., 2023] exhibit irregular geometries represented by unstructured meshes or point clouds, making them incompatible with the regular grids required by spectral operators [Li et al., 2021]; (iii) Long-horizon instability: Autoregressive rollout accumulates errors rapidly, often causing predictions to diverge within 10-20 steps [Brandstetter et al., 2022; Lippe et al., 2023]. These challenges compound: irregular geometries make standard derivative computation unstable, tainting physics residuals and accelerates long-horizon divergence.

Recent research pursues three directions. Spectral neural operators [Li et al., 2021] like FNO achieve efficient global modeling, but depend on regular grids; handling irregular geometries requires coordinate mappings [Li et al., 2023a] that incur significant computational cost. Geometric deep learning [Bronstein et al., 2021] leverages GNNs [Fortunato et al., 2022] and Transformers [Hao et al., 2023] to handle generalized geometries. However, these approaches are impractical when labeled data is scarce or costly. To bypass the dependency on labeled data, Physics-informed learning (PINNs) [Raissi et al., 2019] allows for unsupervised learning using PDE residuals; however, it is prone to memory inflation and gradient pathologies when computing higherorder derivatives [Yu et al., 2022]. Alternatives like variational formulations [Kharazmi et al., 2021] and probabilistic models [Zhang et al., 2025] mitigate oscillations but not the memory cost. However, these methods neglecting temporal structure suffer error accumulation. Crucially, current paradigms impose an trade-off: spectral methods sacrifice geometric flexibility for efficiency, geometric learners demand expensive supervision, and PINNs compromise stability for label-free training. As a result, these requirements remain difficult to jointly achieve within a single model.

To overcome these issues, we explore physics restrictions from a holistic standpoint. Traditional differential residuals require pointwise derivatives [Raissi et al., 2019], which are unstable on irregular geometries. Instead, we reformulate PDE evolution as state propagation via integral consistency. This replaces differential residuals with integral consistency constraints, avoiding derivative computation and enabling label-free learning. This integral formulation, developed for advection-diffusion systems [Karniadakis et al., 2021], enhances stability for parabolic PDEs. We bridge the gap between spectral operators and geometric learning by projecting generalized geometries onto a shared regular latent space, combining spectral efficiency and geometric flexibility.

To realize this vision, we present the Physics-Constrained Unsupervised Neural Operator for Long-Horizon PDE Learning on Generalized Geometries (UNOP). This framework delivers three key contributions:

• Latent Integral Physics Embedding (LIPE): A derivative-free objective enforcing physical laws via stochastic spatiotemporal integral sampling. Designed for advection-diffusion systems, it replaces unstable differential residuals to enable robust unsupervised learning without simulation labels.
• Geometry-Agnostic Latent Adapter (GALA): A unified interface that projects generalized geometries into a regular latent space by aggregating features from grid and their neighboring vertices. This latent representation provides a consistent domain for LIPE, enabling accurate evaluation of spatial integrals while decoupling geometric complexity from operator computation.
• Gated Spectral Evolution Operator (GSEO): A spectral evolution module employing curvature-aware gating on spectral features to preserve spatial regions with high curvature. It accurately captures multiscale dynamics and maintains shocks and discontinuities in the solution field during long-horizon rollouts, preventing the numerical dissipation that affects standard spectral methods.

2 RelatedWork

Current research on neural PDE solvers primarily focuses on three aspects: architecture efficiency, physics constraints, and long-horizon stability. However, achieving all three simultaneously on generalized geometries remains an elusive goal.

2.1 Spectral Methods and Geometric Learning

Spectral methods achieve efficient global modeling. FNO [Li et al., 2021] learns integral kernels with O(N logN) complexity, while DeepONet [Lu et al., 2021] establishes the operator learning framework [Kovachki et al., 2023]. Extensions such as U-FNO [Wen et al., 2022], U-NO [Rahman et al., 2023], CNO [Raonic et al., 2023], and adaptive FNO [Guibas et al., 2021] enhance multiresolution performance but rely on regular grids. Geo-FNO [Li et al., 2023a] and domain-adaptive methods [Lingsch et al., 2024] attempt mapping irregular geometries to regular grids but remain topologically constrained. Conversely, geometric deep learning [Bronstein et al., 2021] handles generalized geometries using GNNs [Sanchez-Gonzalez et al., 2020; Pfaff et al., 2021; Fortunato et al., 2022] or Transformers [Hao et al., 2023; Hao et al., 2024; Wu et al., 2024]. Specialized models like GAOT [Wen et al., 2025] and AeroGTO [Liu et al., 2025] target 3D aerodynamics. While flexible, these methods depend on expensive CFD/FEM simulations for supervision, motivating the search for data-efficient, physics-constrained alternatives.

2.2 Physics-Informed Learning

To minimize supervision, PINNs [Raissi et al., 2019] embed PDE governing equations into loss functions via automatic differentiation, with variational formulations like Deep Ritz [E and Yu, 2018] as alternatives. Subsequent works address gradient pathology [Yu et al., 2022] and propagation failure [Wu et al., 2025]. PINO [Li et al., 2024] combines data and residual losses. However, computing higherorder derivatives on neural outputs amplifies high-frequency noise and causes prohibitive memory costs on 3D geometries [Yu et al., 2022]. Weak-formmethods [Kharazmi et al., 2021] and MCNP [Zhang et al., 2025] offer limited mitigation. MCNP employs Monte Carlo integration but is primarily limited to low-dimensional regular grids due to interpolation constraints. Moreover, these methods often struggle to maintain stability over long time horizons.

2.3 Long-TermPrediction

Error accumulation during autoregressive rollout is a key bottleneck [Lippe et al., 2023]. Existing techniques, including noise injection [Brandstetter et al., 2022], temporal bundling [Brandstetter et al., 2022], and solver-based refinement [Giorgini et al., 2025; Lippe et al., 2023], provide partial mitigation. Pure data-driven methods [Cachay et al., 2023] often drift over long horizons due to the lack of physical constraints. While recent studies incorporate physicsbased regularization [Wang and Perdikaris, 2023], most rely on differential residuals, which are sensitive to the spatial and temporal discretization of the solution, leading to error amplification during long-term rollouts.

2.4 Our Contribution

UNOP unifies these strengths by decoupling geometry from physics. We introduce GALA to project irregular geometries into a regular latent space, enabling efficient spectral processing via GSEO. To ensure stability, we propose LIPE, which replaces unstable differential residuals with integral consistency. By enforcing physical laws via stochastic sampling, LIPE eliminates gradient pathologies, stabilizes predictions, and prevents error accumulation in long-horizon rollouts.

3 Methodology

3.1 ProblemFormulation

We consider physical systems governed by the advectiondiffusion equation [Karniadakis et al., 2021] on a domain Ω ⊂ Rd (d ∈ {1, 2, 3}):

where u(x, t) : Ω×[0, T ] → RC is the state field comprising C physical quantities, v ∈ Rd is the advection velocity, ν ≥

0 is the diffusion coefficient (possibly spatially varying ν(x)), and f is an external source term.

The domain is discretized as N samples {(xi, ui)}Ni=1 , where xi ∈ Rd is the spatial coordinate and ui = u(xi, t) ∈ RC denotes the physical quantities at the i-th point. This discretization may reside on regular grids or irregular geometries, all of which are treated uniformly by our framework.

Our goal is to learn the solution operator G : ut ↦→ ut+∆t in an unsupervised manner without labeled simulation data, maintaining long-horizon stability.

3.2 Framework Overview

As illustrated in Figure 1, UNOP unifies three modules into a single pipeline. Crucially, LIPE (Sec. 3.4) constitutes the core physical principle, while GALA and GSEO provide the numerical realization.

To clarify the data flow and symbol usage, we distinguish three operational spaces:

• Physical Space (x): The original domain containing irregular input points xi.
• Latent Space (Z): A unified representation where generalized inputs are mapped for efficient evolution.
• Integral Space (I): The stochastic domain where LIPE
• constraints are evaluated via interpolation on latent grid. The forward pass is defined as:

where Zt = E(ut) is the encoded latent state, S is the GSEO

evolution operator, and D is the decoder. The pipeline is trained end-to-end via the LIPE objective to solve Eq. (1).

3.3 GALA: Geometry-Agnostic Latent Adapter

Standard spectral methods require regular grids, but realworld data often comes as irregular point clouds. GALA projects the discretized PDE data {(xi, ui)}Ni=1 onto a unified regular latent grid, enabling efficient spectral processing and stable interpolation for LIPE.

Local Feature Extraction

To accommodate the discretized PDE data {(xi, ui)}Ni=1 shown in Figure 1, we employ an encoding strategy. We extract local features hi ∈ RD based on the spatial arrangement of the points:

For regular inputs, we use a lightweight CNN [Krizhevsky et al., 2012] to extract local spatial features, while EGraph is implemented as a physics-token based GATv2 [Brody et al., 2022] to handle complex topologies. For the graph branch, ui and xi denote the physical quantities and coordinates of point i, aggregating information from neighbors j ∈ Ni. This step ensures that local geometric features are extracted before unified projection.

Latent Grid Projection

To construct a structured representation, features {hi} are projected onto a regular latent grid G with Ng vertices

{gk}Ng k=1 . Each vertex aggregates nearby input points in the neighborhood Nk, weighted by a learnable kernel ϕ based on spatial distance:

Here, Z(gk) is the feature vector at vertex gk, corresponding to a component of the latent state Zt. The resulting tensor Zt serves as the uniform substrate for GSEO evolution. The kernel ϕ is based on a Gaussian approximation with learnable parameters (see Appendix for implementation and weight derivation).

Importantly, the dimensionality of this latent grid does not correspond to the physical dimensionality of the input domain. Regardless of whether the original PDE is defined in 1D, 2D, or 3D, GALA embeds the geometry into a unified latent representation Zt, enabling a shared evolution mechanism. Furthermore, the number of grid vertices Ng and the latent feature dimension D are hyperparameters independent of the input size, allowing the framework to scale to highresolution problems without architectural changes.

3.4 LIPE: Latent Integral Physics Embedding

LIPE defines the unsupervised training objective. Instead of minimizing unstable differential residuals, we employ integral consistency derived from the physics of diffusion, enforced directly on the decoded prediction ˆu = D(Z) as shown in Eq. (2).

Integral Targets

Drawing inspiration from the Feynman-Kac formula [Karniadakis et al., 2021] which links PDEs to stochastic processes, we reformulate Eq. 1 into integral forms.

Temporal Target: For transient dynamics, the state at t + ∆t is theoretically given by the expectation:

Here, x - v∆t represents deterministic advection, ξ ∼ N (0, 2ν∆t Id) represents Brownian motion diffusion [Zhang et al., 2025], and f(x)∆t accounts for the integration of the external source term f in Eq. (1).

Spatial Target: For steady-state (∂tu = 0), mass conservation dictates zero net flux across any control volume Ωr:

where u is the velocity vector, n is the outward unit normal vector on the boundary ∂Ωr, and dS denotes the surface area element.

Stochastic Discretization

To compute these targets during training, we define the integral operator I via stochastic sampling.

For Temporal Dynamics: We approximate the expectation Eξ in Eq. 5 using K stochastic samples:

where ξ(k) are sampled noise vectors and wk are Gaussian weights normalized to sum to unity. Crucially, evaluating u at the off-grid locations x - v∆t + ξ(k) is enabled by interpolation on the GALA latent grid.

For Spatial Constraints: We approximate the surface integral via discrete summation:

where sj are sampled points on the spherical surface Sr, nj are the corresponding outward normals, r is the radius of Sr, and M is the number of surface samples. Here, ˆut+∆t = D(Zt + ∆t S(Zt)) denotes the decoded prediction for the next time step. In practice, we select K and M to balance computational cost and variance reduction, with ablation studies in Sec. 4 confirming stability.

Training Objective

While the integral targets are derived from physical variables u, the LIPE objective is computed on the decoded network prediction ˆu = D(Z). This formulation ensures that the supervision signal from physical laws backpropagates through the decoder to guide the latent evolution in GALA and GSEO. The temporal loss is:

where the term f∆t represents the explicit integration of the external source field f, consistent with Eq. 5.

For spatial problems, we minimize the flux residual:

This loss term explicitly penalizes local mass violations. Crucially, these losses provide the only supervision signal, guiding the entire framework to discover physical dynamics from scratch.

3.5 GSEO: Implementing Latent Evolution

GSEO provides a differentiable evolution prior in latent space. Without an explicit evolution operator, LIPE alone would lead to ill-posed optimization, as infinitely many mappings can satisfy integral constraints locally while violating global dynamics.

Spectral Processing

Global mixing is performed in the frequency domain:

where F is the Fourier transform and R are learnable weights. This effectively captures global advection features.

Gated Update with Curvature Awareness

Latent states corresponding to generalized geometries often contain high curvature regions, such as shocks or sharp gradients, which are common in realistic physical systems. To address this, we introduce a curvature-aware gate Γ that modulates the spectral evolution:

x

x

x

x

x

Figure 2: Qualitative Results on Spatiotemporal Accuracy. Left: Absolute error evolution for 1D Convection-Diffusion. UNOP maintains negligible error levels. In contrast, PINO and PI-DeepONet exhibit significant deviations; their error curves are frequently truncated at the visualization limit, appearing as saturated spikes that far exceed the unified scale. Right: Error maps for 2D Navier-Stokes. UNOP exhibits a uniformly deep blue distribution indicating low error. Conversely, SP-FNO shows distinct striping artifacts caused by spectral leakage, and MCNP displays red high-error regions in turbulent zones where it fails to resolve fine scales.

The gate Γ is defined as a function of both geometric and physical curvature indicators:

Here, ¯κgeom(Zt) characterizes geometric curvature induced by the spatial arrangement of input points, while ¯κphys(Zt) is the Laplacian magnitude of the latent state, reflecting sharp physical variations. The learnable parameters α, β, and γ balance their contributions, enabling adaptive modulation of the update strength in high-curvature regions and mitigating the smoothing artifacts of standard spectral methods.

3.6 Decoding

The final prediction is obtained by decoding the evolved latent state Zt+∆t:ˆut+∆t = D(Zt+∆t) + ψ(h) (14)

Here, the decoder D reconstructs physical values from the latent grid using interpolation, and may incorporate a residual network that leverages the pre-projection features h (Eq. 3) to recover high-frequency details lost during grid projection.

In summary, each component of UNOP addresses a specific limitation: GALA handles generalized geometries, LIPE enforces stable and consistent evolution over time, and GSEO for preserving fine-scale details in long-horizon predictions. The effectiveness of this design is confirmed in Sec. 4.

4 Experiments

Datasets and Tasks. We evaluate UNOP on three benchmarks, each specifically chosen to validate a core component of our framework: (1) 1D Convection-Diffusion: This benchmark tests UNOP's long-horizon stability on abrupt changes, highlighting LIPE's role in managing discontinuities. We define six regimes based on advection velocity β and diffusion coefficient κ: Regimes E1-E3 use low advection β = 0.01 with diffusion κ ∈ {0.001, 0.005, 0.01} respectively. Regimes E4-E6 use high advection β = 0.1 with the same diffusion sequence.(2) 2D Navier-Stokes: This chaotic system evaluates GSEO's ability to capture complex dynamics and maintain stable long-term predictions. Regimes E1 and E2 use Zero forcing, E3 and E4 use Li forcing, and E5 and E6 use Kolmogorov forcing. Within each forcing type, the odd-numbered regime uses high viscosity ν = 10-4 , while the even-numbered regime uses low viscosity ν = 10-5 .(3) 3D ShapeNet-Car: A spatial domain task on unstructured point clouds, aimed at demonstrating GALA's adaptation to generalized geometries.

Comparison Methods. UNOP is compared with state-ofthe-art unsupervised neural operators: PI-DeepONet [Wang and Perdikaris, 2023] serves as a physics-informed baseline for 1D functional mapping. PINO [Li et al., 2024] combines data and residual losses within an FNO architecture. SP-FNO [Tran et al., 2023] utilizes a factorized FNO variant for efficient 2D modeling. MCNP [Zhang et al., 2025] acts as a integral-based competitor using Monte Carlo integration. Implementation Details. All models are implemented in PyTorch and trained on single NVIDIA 4090 GPU. Detailed hyperparameters are provided in Appendix. All quantitative results are reported as the mean of three independent experiments.

Evaluation Metrics. We assess performance using Relative L2 Error. For 3D Car, we report the Drag Coefficient Error Cd and physics violation metrics EBC and EMass.

4.1 Main Results

Results on 1D/2D Benchmarks. Tables 1 and 2 present the quantitative comparison. In Table 1, regarding the 1D Convection-Diffusion benchmark, UNOP achieves the lowest error across all regimes. Notably, in the high-advection regime E4 (β = 0.1, κ = 10-3 ), UNOP reduces the error to 0.09%, outperforming PINO at 3.32% and PI-DeepONet at 8.60% by orders of magnitude. This improvement mainly comes from LIPE's integral formulation, which successfully eliminates the high-frequency oscillations inherent in differential baselines. Table 2 details the results for the 2D NavierStokes benchmark, where UNOP consistently outperforms all other baselines. For instance, in the Kolmogorov forcing task E6, UNOP achieves 3.02% error at rollout 20, significantly lower than SP-FNO at 13.6% and MCNP at 9.95%. This confirms that GSEO's curvature-adaptive gating effectively mitigates spectral leakage, preserving fine-scale turbulent structures that pure spectral methods tend to smooth out.

Spatiotemporal Fidelity. Figure 2 provides a qualitative comparison. In 1D tasks, UNOP achieves near-zero deviation from the ground truth shock, with mean squared error around 10-4 , while baseline methods exhibit noticeable spurious oscillations. In 2D tasks, UNOP maintains low error across the entire domain, accurately capturing smooth low-frequency patterns while preserving high-frequency structures such as sharp gradients and vortices. In contrast, SP-FNO produces striping artifacts, and MCNP fails in high-gradient regions.

Long-term Stability Analysis. Figure 3 and Figure 4 analyze error evolution during extended rollouts. As shown in Figure 4, during 1D extrapolation at T = 50, baseline errors surge: PI-DeepONet diverges beyond 120% relative error, while MCNP drifts up to 20%. In contrast, UNOP maintains a stable error trajectory below 15% throughout the rollout. LIPE stabilizes long-term rollouts by enforcing integral consistency, effectively suppressing autoregressive error accumulation. In 2D turbulence, UNOP preserves sharp vortex structures at T = 20 in Figure 5, whereas MCNP blurs fine details. This performance arises from a careful balance of global and local dynamics: low-frequency convective patterns are retained to stabilize the rollout, while GSEO's curvature adaptive gating selectively injects high frequency residuals into regions with sharp gradients. By doing so, the method recovers fine scale details lost to numerical dissipation in purely spectral or interpolation based approaches.

0

Table 3: 3D Aerodynamics on ShapeNet-Car. Comparison under Semi-Supervised (10% Label) settings. EBC/Mass: Physics violations (×10-3 ). L2: Relative error (%).

Results on 3D Generalization. Table 3 reports the performance on the semi-supervised ShapeNet-Car task using 10% labels. UNOP achieves an L2 error of 18.97% and a Cd error of 4.06%, significantly outperforming the best baseline MCNP by 30%. The methods such as PINO and PIDeepONet exhibit substantial performance degradation with errors exceeding 47%. This highlights the inherent limitation of standard operator architectures when applied to irregular geometries without dedicated geometric adaptation, whereas UNOP's success validates GALA's ability to decouple geometric complexity from physical learning. UNOP with 10% supervision outperforms the fully supervised UNOP shown in gray, suggesting that sparse supervision combined with integral physics regularizes overfitting.

4.2 Ablation Studies

Table 4 validates the contribution of each component. Effect of LIPE. Replacing LIPE with standard differential loss increases the 2D error from 3.25% to 33.13%. This dramatic 10× error increase highlights LIPE's critical role in UNOP's stability. Effect of GALA. Removing GALA in 3D tasks leads to training failure with errors exceeding 100%, demonstrating its necessity for handling irregular geometries. Effect of GSEO. Disabling the Curvature Gate increases the error to 23.09% on 2D tasks, indicating its importance in preserving sharp dynamics that spectral filters tend to smooth.

4.3 3DVisualization Analysis

3D Boundary Adherence. Figure 6 visualizes surface pressure and velocity. UNOP maintains strict boundary adherence, whereas baselines show significant violations near boundaries, further confirming the efficacy of our integral boundary constraints.

5 Conclusion

We propose UNOP, a framework for unsupervised PDE learning that enforces integral consistency via LIPE rather than differential residuals. By reformulating physical constraints in latent space, UNOP achieves a stable balance between geometric flexibility and temporal evolution, demonstrating superior accuracy and long-term stability across 1D-3D benchmarks while capturing both global patterns and local structures such as shocks and vortices. Via the GALA interface, it handles complex geometries in high dimensions, consistently outperforming interpolation and spectral baselines in preserving physical fidelity over extended rollouts.

Limitations and Future Work. The current LIPE formulation targets advection-diffusion PDEs via the Feynman-Kac framework. Extending it to hyperbolic systems with strong shocks or multi-physics couplings remains an important direction for future work.

Contribution Statement

Xinrui Cheng and Ruihan Hu contributed equally to this work. ∗Ruihan Hu is the corresponding author.

Ethical Statement

There are no ethical issues.

Acknowledgments

This work was supported in part by the National Key Research and Development Program under Grant 2023YFF0905500, in part by the National Key Research and Development Program of Heilongjiang under Grant 2022ZX01A11, in part by the Natural Science Foundation of Heilongjiang Province under Grant LH2023F001, in part by the ZTE Industry-University-Institute Cooperation Funds under Grant No. 2022H2ZTE02-01-P4, and in part by the Spring Goose Talent Team Program of Heilongjiang.