arXiv CSopen access

GAIA: Geometry-Adaptive Operator Learning for Forward and Inverse Problems

GAIA: Geometry-Adaptive Operator Learning for Forward and Inverse Problems∗ Meenakshi Krishnan1 , Pranav Pulijala1† , Ke Chen2 , Haizhao Yang1 , and Ramani Duraiswami1

arXiv:2607.01128v1 [cs.LG] 1 Jul 2026

1

University of Maryland, College Park 2 University of Delaware, Newark

Abstract Operator learning for partial differential equations (PDEs) on arbitrary geometries builds fast neural surrogates for large-scale simulation. Although recent geometryadaptive neural operators have made substantial progress, they are mainly designed for forward problems in which inputs and outputs share the same spatial domain. This limits their applicability for boundary value problems (BVPs) and inverse problems, where inputs and outputs may live on different domains. We introduce the Geometry-Adaptive Integral Autoencoder (GAIA), an operator learning model that encodes the domain boundary and the interior field distribution into geometry tokens, and conditions integral transform layers on these tokens via cross-attention, allowing the kernel to adapt locally to geometric features. This yields a single architecture for forward (including BVPs) and inverse problems on arbitrary domains in one pass, without retraining, iterative optimization, or graph construction. We evaluate GAIA on seven 2D and 3D benchmarks, four of which are new or substantially extended benchmarks for inverse problems and BVP: electrical impedance tomography, optical tomography, 3D Darcy flow on varying geometries, and a modified setting of Poisson BVP on mechanical components benchmark (MCB). GAIA sets new state-of-the-art results on every inverse and BVP task, reducing median relative L2 error by 64% on airfoil flow reconstruction and 27% on EIT relative to the next best amortized method, and outperforming all baselines on every shape category of MCB. On other forward problems, GAIA is competitive with specialized solvers while maintaining stable accuracy across point resolutions on which transformer-based baselines degrade.

1

Introduction

Partial differential equations are used to model problems across science and engineering, from fluid dynamics [1] to medical imaging [2]. Traditional solvers are accurate but expensive, particularly in multi-query settings such as inverse problems, optimal design, and uncertainty quantification, where the forward model is invoked thousands of times. Early neural methods developed to address these costs, such as Physics-Informed Neural Networks (PINNs) [3] and deep inverse priors [4] are instance-specific, requiring retraining for every new boundary condition or parameter change. To overcome the need for repeated retraining, operator learning models instead learn maps between function spaces directly. However, classical neural operators such as FNO [5] and DeepONet ∗ This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department

of Energy Office of Science User Facility, operated under Contract No. DE-AC02-05CH11231 using NERSC award ALCC-ERCAP0034775. † Significant technical contribution. Correspondence: Meenakshi Krishnan ([email protected]). Preprint.

[6] require regular grids or fixed discretizations. Many problems in structural analysis, medical tomography, and aerodynamic design are posed on domains whose geometry varies across problem instances, in both forward and inverse settings. A broadly applicable neural proxy must thus apply to arbitrary varying geometries without retraining. Recent geometry-adaptive operators, such as GAOT [7], GINO [8], Transolver [9], handle unstructured point clouds via graph neural networks or transformer attention. However, these methods are mostly designed for forward modeling. Some, such as Transolver, require the input and output to share the same spatial discretization, precluding problems where the input and output live on different point sets - as in BVPs and inverse problems with boundary measurements. Applying these architectures to the latter requires computationally expensive iterative optimization loops over forward solvers. Others, such as GAOT, are designed and evaluated exclusively for forward problems. While GINO can in principle handle inverse problems, its reliance on graph message passing incurs substantial memory and computational overhead. We introduce the Geometry-Adaptive Integral Autoencoder (GAIA), a unified model for both forward and inverse problems on arbitrary domains. Building on IAE-Net [10], we introduce a dual-pathway tokenization that encodes both the domain boundary and the interior field distribution, and condition the integral kernels on these tokens via multi-head cross-attention. Unlike methods that encode geometry implicitly through per-sample neural field fitting [11] or require expensive graph construction [8, 7], this provides explicit, spatially adaptive geometric context at each integration point. The query structure accommodates problems where the input and output live on different spatial domains, enabling efficient single-pass solution of inverse problems without iterative optimization. In summary: • Unified forward/inverse geometry-adaptive model. We propose GAIA, a geometryadaptive integral operator that conditions learned kernels on domain geometry via crossattention, enabling single-pass solution of forward and inverse problems on arbitrary geometries. • New benchmarks for varying-geometry inverse and BVP problems. We introduce new varying-geometry benchmarks including for inverse problems such as Electrical Impedance Tomography (EIT), Optical Tomography (OT); all of which will be made publicly available. • Empirical validation. We validate GAIA on seven time-independent PDE benchmarks spanning 2D and 3D problems. GAIA achieves state-of-the-art on all inverse and BVP tasks, reducing median relative L2 error by 27% on EIT and 64% on airfoil relative to the next best method, while remaining competitive on forward benchmarks with stable accuracy across resolutions where transformer baselines degrade.

2

Related Work

2.1

Operator Learning

Early scientific machine learning methods such as Physics-Informed Neural Networks [3, 12–16] and CNN-based surrogates [17–19] are either instance-specific or restricted to regular grids; we focus on operator learning. Neural operators learn mappings between function spaces. DeepONet [6] uses a branch-and-trunk decomposition but requires inputs sampled on a fixed grid. FNO [5] parameterizes integral kernels in the spectral domain via the Fast Fourier Transform (FFT), achieving good performance on regular grids. Numerous extensions improve scalability and discretization handling through low-rank or factorized kernels [20–23], graph message passing [24, 25], and wavelet bases [26]. IAE-Net [10] achieves discretization invariance through a stack of integral autoencoder blocks with data-driven kernels, but its kernels depend only on local coordinates and function values with no access to global domain geometry. For inverse problems, neural operators have been used either as differentiable surrogates within iterative optimization, MCMC, or diffusion-based sampling loops [27–31], or to directly learn amortized inverse maps from measurement-to-parameter data pairs, enabling single-pass inference [32, 33]. The former does not avoid the cost of iterative sampling at inference. GAIA belongs in the latter group. Neural Inverse Operators [32] learn operator-to-function mappings via compositions of DeepONet and FNO; invertible architectures such as iFNO [33] jointly learn forward and inverse 2

Boundary

Boundary Points

Tokenizer

Mesh Nodes Field Values

State Slice Tokenizer

Boundary Tokens

LayerNorm Geometry + Tokens GeLU

⊕

Slice Tokens

Transformer Encoder (Token Updaters)

Output Field Field Values Input

FC Features

Adaptive Block 1

Adaptive IAE Block 2

· …

Adaptive IAE Block K

Reducer (Conv1D)

FC

Adaptive Encoder

Norm & FeedFwd

Integral Transform (x^T @ kernel)

Expand

Mesh Nodes

⊕

Merge (Add)

Output

Linear

Norm & FeedFwd

Integral Transform (x^T @ kernel)

MLP

Softmax

Geometry Tokens

Token Attention

NeRF MultiHead Positional Attention MLP Encoding

Mesh Nodes

⊕

Linear

Geometry Tokens

NeRF MultiHead Positional Attention MLP Encoding

Token Attention

Features

Compressor (Linear)

AdaptiveIAEBlock

Adaptive Decoder

Skip Connection (Conv1D)

Figure 1: GAIA architecture. Top: Two complementary tokenizers produce geometry tokens from boundary points and interior mesh-field values, refined by transformer encoder layers. Middle: A stack of geometry-conditioned Adaptive IAE blocks with DenseNet-style skip connections maps input features to the output field. Bottom: Each block follows an encode–process–decode pattern with integral-transform kernels conditioned on geometry tokens via cross-attention. Separate encoder and decoder query sets enable inverse problems and BVPs.

maps through shared invertible blocks. However, both are tied to fixed regular grids due to their reliance on FFT. See [34] for a probabilistic perspective on operator learning for inverse problems. 2.2

Geometry Adaptive Operators

Geometry-adaptive operators take three broad approaches. Coordinate-transform methods like Geo-FNO [35] learn a deformation from the physical domain to a regular latent grid; this works for topologically simple domains but fails on complex shapes. Graph- and message-passing methods including GINO [8], GAOT [7], and RIGNO [36] handle unstructured point clouds natively but require per-sample graph construction that dominates inference cost. Implicit-representation methods [37, 38] include CORAL [11], which encodes fields as implicit neural representations and learns maps between latent codes. Transolver [9] and LNO [39] are attention-based methods which project mesh points into a small set of learnable physical states via attention; broader transformer operators include GNOT [40] and OFormer [41, 42]. There are also specialized architectures for specific PDE such as Neural Greens Function [43] for symmetric linear PDEs. Among these, only LNO has been demonstrated on inverse tasks with differing input-output domains; CORAL and GeoFNO demonstrate iterative optimization over their forward surrogates for inversion. GAIA differs from all three: unlike CORAL it does not require per-sample neural-field fitting; unlike GINO it constructs no graph; and unlike LNO, which uses attention to compress spatial information into a fixed set of latent states, GAIA conditions the integral kernel itself on geometry at each query point via cross-attention. Our slice tokenizer adapts the soft-clustering aggregation of Transolver [9]; the boundary tokenizer and the use of both token pathways to condition integral kernels via cross-attention are novel, and the results demonstrate the architecture’s efficacy.

3

Methods

We formalize operator-learning, recap the IAE-Net model (§3.1), and describe GAIA (§3.2). 3

3.1

Problem Formulation

We consider a generic time-independent PDE, D(c, u) = q,

∀x ∈ D ⊂ Rd ,

B(u) = ub ,

x ∈ ∂D,

(1)

where u : D → Rm is the solution, c ∈ C represents the physical parameters, q ∈ F is the forcing term, ub prescribes the boundary values; and D and B are the differential and boundary operators, respectively. Also, let χD denote the domain indicator. Here U, C, and F denote suitable function spaces on D. Encapsulating the problem setup into a = (c, q, ub , χD ) ∈ A, the forward operator S : A → U maps a 7→ u. In the inverse setting, we do not have access to the full solution field u. We only have access to observables, often measured at the domain boundaries. This is formalized by a measurement operator M : U → V mapping u to an observation v = M(u). The generalized forward operator F = M ◦ S : A → V maps parameters directly to observables, and the inverse operator F † : V → C recovers c = F † (v, χD ) from these measurements. Throughout, let Ψ : X → Y denote the target mapping (Ψ ≡ S in the forward case, Ψ ≡ F † in the inverse case). Our goal is to approximate Ψ given access to data pairs (f (i) , g (i) ) with g (i) = Ψ(f (i) ). 3.2

IAE-Net

IAE-Net [10] learns an operator Ψ : X → Y between function spaces on compact domains Ωx ⊂ Rdx and Ωy ⊂ Rdy , from finite discretizations Sx ⊂ Ωx , Sy ⊂ Ωy whose cardinalities may vary across samples. The architecture is a recursive stack of L densely-connected autoencoder blocks, performing a series of transformations f → a0 → a1 → · · · → aL → g, with each intermediate function ai defined on Ωa = [0, 1]da . Each block performs an encode–process–decode transformation a → v → u → b, where v and u live on a fixed latent domain Ωz = [0, 1]d , discretized on a grid Sz = {zj }m j=1 independent of the input: Z v(z) =

Z ϕ1 (a(x), x, z; θ1 ) a(x) dx, u(z) = ϕ0 (v(z); θ0 ), b(y) =

Ωx

ϕ2 (u(z), y, z; θ2 ) u(z) dz. (2) Ωz

Here, ϕ1 is the encoding kernel, ϕ0 a latent-space MLP, ϕ2 the decoding kernel, and y ∈ Ωa . The integrals are approximated as discrete sums, mapping variable-cardinality inputs to the fixed grid Sz . This allows IAE-Net to natively handle different input resolutions. Blocks are connected with DenseNet-style skip connections. The kernels ϕ1 , ϕ2 in equation 2 of IAE-Net depend only on local coordinates and the function values a(x), with no access to global domain geometry. For elliptic PDEs, the Green’s function governing the solution operator depends on the domain shape, so a geometry-agnostic kernel cannot capture this dependence. As a result, the same kernel is applied regardless of the domain geometry. Additionally, while IAE-Net’s integral transform structure can in principle support differing input and output domains, this has not been explored for boundary-to-interior mappings. 3.3

Geometry-Adaptive Integral Autoencoder (GAIA)

GAIA conditions IAE-Net’s integral kernels on the domain geometry. We describe tokenizing the boundary points and interior fields, the conditioning mechanism, and the resulting blocks for forward and inverse problems. To address the geometry-agnostic limitation of IAE-Net described in §3.1., the kernels are conditioned on a tokenized representation of the actual domain boundary and interior field distribution at inference time, with each spatial query point attending selectively via cross-attention, providing geometric context that is spatially adaptive. Decoupling query points from integration points lets the encoder integrate over Xbnd while the decoder projects onto Xdom , enabling single-pass solve of inverse problems and BVPs. 3.3.1

Geometry Tokenization

The geometry conditioning is constructed from two complementary token pathways: boundary feature tokens, which encode the domain shape from boundary coordinates alone, and state slice tokens, which encode the joint distribution of spatial coordinates and field values in the domain. Let 4

Boundary Tokenizer

Points

Linear

Layer Norm

Global Max Pool

ReLU

Linear

Layer Norm

Linear

ReLU

Reshape

Tokens

Projection MLP

PointNet Boundary

Nodes

Linear

Softmax

MatMul

Normalize

Linear

Tokens

GeLU

Weights

Linear

Weights

Mesh

Aggregation

SlicerSlicer Slicer

Layer Norm

Features

Embedding ⊕

Field Values

(B, N, H)

2x

2x

Slice Tokenizer

Figure 2: The two tokenizer pathways: GAIA encodes geometry through two complementary token sets. Top – Boundary tokenizer: Boundary point coordinates are processed by a PointNet-style encoder with shared MLPs and global max-pooling permutation-invariant tokens that summarize the global domain shape. Bottom – Slice tokenizer: The interior mesh coordinates and their associated field values are jointly embedded and assigned to soft clusters; mass-normalized aggregation produces tokens that summarize the spatially varying domain physics. Xbnd ∈ RM ×d be the boundary/sensor coordinates and Xdom ∈ RN ×d the interior mesh; let U be the associated input field or measurements. The two pathways produce tokens that are concatenated into a unified T ∈ RK×dt ; K is the total token count, dt the token dimension. Boundary feature tokenization. The boundary coordinates Xbnd are mapped to Kbnd tokens Tbnd ∈ RKbnd ×dt using a PointNet-style encoder [44], which applies shared MLPs followed by permutationinvariant global max-pooling, then projects and reshapes the resulting vector into Kbnd tokens:   Tbnd = Proj max MLP(p) . (3) p∈Xbnd

State slice tokenization. Adapting the soft-clustering aggregation of Transolver [9] to inputs that include both spatial coordinates and physical values, we let X ∈ RP ×d be a point cloud (instantiating as Xdom for forward problems or Xbnd for inverse problems) and U the associated values. The concatenated state [X, U ] is mapped to Kslice soft-cluster tokens via: W = Softmax(MLPcluster ([X, U ])), H = MLPfeat ([X, U ]), (4) Tslice = diag(W ⊤ 1)−1 W ⊤ H, (5) with mass-normalized aggregation; 1 ∈ RP is a vector of ones. A linear projection brings Tslice to the standard token dimension dt . The initial token set T (0) = [Tbnd , Tslice ] is normalized and passed through GELU; between Adaptive IAE blocks, tokens are refined by standard transformer encoder layers, so that each block sees an updated token set T (l) . 3.3.2

Geometry-conditioned integral transforms

GAIA modifies the encoding and decoding integral transforms (Eq. 2) so that the learnable kernel is conditioned on the geometry tokens T (l) . Spatial coordinates are first lifted via Fourier feature encodings γ(x) ∈ R2dF +d with F frequency bands [45]. To inject geometric context into the kernel, we define a spatially-queried context vector C(x) via cross-attention between the encoded coordinates and the geometry tokens:   (γ(x)WQ )(T (l) WK )⊤ √ C(x) = Softmax (T (l) WV ), (6) dk where WQ , WK , WV are learnable projections and dk the scaling factor. The kernel is parameterized by an MLP taking the concatenated spatial features and geometric context: K(x, z; θ) = MLPkernel ([γ(x), γ(z), U (x), C(x)]) , (7) where U (x) is the input feature at point x. The integrals in (2) with this kernel are approximated as: 1 X I(z) = K(xj , z; θ) U (xj ), (8) |X| xj ∈X

followed by a learned residual connection and GELU non-linearity. 5

3.3.3

Adaptive IAE Block

GAIA uses these geometry-conditioned integral transforms as building blocks. Each Adaptive IAE block follows the encode–process–decode structure (equation 2) of IAE-Net, but with kernels conditioned on geometry tokens via cross-attention (Equation 7). The encoder maps from the N -point input to m fixed latent modes, a pointwise MLP processes the latent representation, and the decoder maps back to N points using the same conditioned kernel queried at the output coordinates. This bottleneck gives per-block complexity O(N m) with m ≪ N . 3.3.4

Main architecture

For forward problems where input and output share the same spatial domain, the Adaptive IAE blocks perform volume integrals over the domain mesh. For inverse problems and BVPs, an additional module first transfers information from boundary measurements to the interior. Forward problems. The input field U ∈ RN ×c is defined on Xdom . The network applies L Adaptive (l) IAE blocks with dense skip connections: the input Hin to the l-th block is the concatenation of all preceding outputs, compressed back to width w via 1D convolution. The final prediction comes from passing the densely concatenated outputs through a reduction layer and a GELU MLP. Inverse problems and BVPs. Here the input is sparse measurements U ∈ RM ×cobs at boundary locations Xbnd , and the target lives on the dense interior mesh Xdom . To bridge this gap, an Observationto-Domain Decoder first transfers information from the boundary to a fixed set of Kmodes learnable latent coordinates Z ∈ RKmodes ×d , then projects from the latent set onto the interior mesh. Concretely, the boundary data is embedded as Hbnd = MLPin (U ) ∈ RM ×w . A geometry-conditioned boundary integral lifts it to the latent set: 1 X Hmode (zk ) = Kenc (xj , zk ; T (0) ) Hbnd (xj ) ∀zk ∈ Z. (9) |X| xj ∈Xbnd

A fixed-size MLP processes the latent representation, H̃mode = MLP(Hmode ), after which a second geometry-conditioned integral projects onto Xdom : 1 X H (0) (yi ) = Kdec (yi , zk ; T (1) ) H̃mode (zk ) ∆zk , ∀yi ∈ Xdom . (10) |Z| zk ∈Z

The resulting H (0) ∈ RN ×w enters the stack of Adaptive IAE blocks identically to the forward case.

4

Results

We evaluate GAIA against five geometry-adaptive baselines—GINO [8], CORAL [11], Transolver [9], GAOT [7], and LNO [39]—on seven benchmarks spanning forward, inverse, and BVP settings. Four of these benchmarks are new or extended contributions (Section 4.1); the remaining three are standard benchmarks. Transolver requires shared input/output discretization and is omitted from inverse and BVP comparisons. NIO [32], the closest amortized inverse-operator method, is FNO-based and operates on fixed regular grids; it cannot be applied to our varying-geometry benchmarks by construction. For the 3D Poisson BVP we additionally compare against NGF [43], which is restricted to forward problems for linear symmetric PDEs. All models are trained on a single NVIDIA A6000 GPU. We minimize the relative L2 error,   ∥(Ψn (f ; θ) − Ψ(f )) ⊙ M ∥2 L(θ) = Ef ∼µ , (11) ∥Ψ(f ) ⊙ M ∥2 + ϵ where ⊙ denotes the Hadamard product. We use masking for varying mesh sizes; M is a binary mask indicating valid interior nodes. Optimizers, learning rates, schedules, and full architectural hyperparameters for both GAIA and baselines are reported in Appendix B. 4.1 4.1.1

Benchmarks Inverse problems

Electrical Impedance Tomography (EIT). EIT is an imaging technique that reconstructs the internal conductivity a(x) from boundary electrical measurements, solving an ill-posed [46] inverse problem 6

Figure 3: Model predictions on inverse benchmarks. Top: EIT conductivity reconstruction. Bottom: Airfoil Mach field reconstruction from sparse noisy observations. GAIA captures fine-scale features like the conductivity patterns and transonic shock that baselines either mislocate (GINO) or diffuse (CORAL, LNO).

governed by ∇ · (a(x)∇u(x)) = 0 on Ω ⊆ R2 . We apply L=20 Dirichlet excitation patterns at M =272 sensors on star-shaped domains and measure the Neumann response; the inverse operator maps the resulting Dirichlet-to-Neumann data Λa to a(x). We extend the fixed-geometry setting of [32] to varying star-shaped domains. Optical Tomography (OT). OT is an imaging technique that recovers the scattering coefficient σs (x) of a tissue-like medium from boundary light measurements. The transport of photons through the medium Ω ⊂ R2 is governed by the stationary radiative transport equation (RTE): Z Φ(v · v ′ ) u(x, v ′ ) dv ′ , x ∈ Ω, (12) v · ∇x u(x, v) + σt (x) u(x, v) = σs (x) S d−1

where u(x, v) is the photon density at position x traveling in direction v, σt = σs + σa is the total attenuation, and Φ is the scattering phase function. The boundary observable is the Albedo 1 1 operator Λ R : L (∂Ω) → L (∂Ω), which maps incoming illumination ϕ(x) to the outgoing flux J+ (x) = v·nx >0 (v · nx ) u(x, v) dv. The forward-peaked anisotropic scattering is modeled via the Henyey-Greenstein phase function. The Albedo operator is discretized as a dense source-to-receiver matrix. We use difference imaging [47]: the Albedo operator Λ0 of a homogeneous baseline (σs,0 =1) is used as a reference, and the inverse operator maps the difference ∆Λ = Λ − Λ0 to the interior scattering perturbation δσs (x). Domains are random convex pentagons. This is a new dataset. Airfoil reconstruction. We use the transonic airfoil dataset of [35], which contains steady-state solutions of the compressible Euler equations at M∞ = 0.8 over randomly deformed NACA-0012 profiles. We reformulate a sparse-to-full inverse reconstruction task: only 10% of mesh nodes are observed, and at inference, the observations are corrupted by 1% relative Gaussian noise. The model uses the sparse noisy observations to recover the full Mach field over the domain. 4.1.2

Forward problems

3D Poisson BVP on mechanical components. We solve −∆u = f on 3D mechanical parts from the MCB dataset [48], spanning four shape categories: gears, nuts, fittings, and screws & bolts. We adapt the setting and meshes of [43], but formulate a pure BVP with fixed source and varying randomized Dirichlet boundary conditions, solved with FEniCS [49]. The target operator thus maps the Dirichlet boundary function to the solution on the interior. 3D Darcy flow. We learn the map from a log-normal permeability field κ(x) to the pressure u(x) satisfying −∇ · (κ∇u) = 1 with homogeneous Dirichlet conditions on star-shaped 3D domains defined by random spherical harmonic expansions. Poisson-Gauss and Elasticity. These are standard benchmarks with shared input–output discretizations from [36] and [35] respectively. Details are in Appendix A. 7

Table 1: Comparison on inverse problems. GAIA has the best performance across all benchmarks. Median relative L2 error [%] (point cloud) Dataset

GINO

CORAL

LNO

GAOT

GAIA

EIT OT Airfoil

28.73 1.79 4.50

1.66 2.93 1.60

0.97 1.82 2.26

7.69 5.91 9.95

0.71 1.41 0.58

Table 2: BVP on Mechanical Components Benchmark. Median L2 relative error. GAIA achieves best error across all shapes.

4.2 4.2.1

Model

Fitting

Gear

Nut

Screws/Bolts

CORAL NGF LNO GAIA

18.82 11.66 14.25 10.15

13.57 5.89 6.14 3.44

23.25 5.48 6.50 5.16

20.76 6.52 5.52 3.68

Discussion Inverse Problems

Table 1 reports median relative L2 error on inverse and BVP benchmarks; GAIA achieves the lowest error on every task. Against the next best amortized method, GAIA reduces error by 27% on EIT (0.71% vs. 0.97% for LNO), 21% on OT (1.41% vs. 1.79% for GINO), and 64% on the Airfoil task (0.58% vs. 1.60% for CORAL). The Airfoil margin is the largest, likely because baselines struggle with the combination of noise, sparsity, and geometric variation. Visual reconstructions in the Appendix (Figures 7-13) show GAIA recovering fine-scale features (such as oscillatory conductivity patches in EIT, transonic shock structure in Airfoil) that baselines either smear or miss entirely. Table 3: Comparison on forward problems. GAIA is competitive with SoTA methods. Median relative L2 error [%] (point cloud)

4.2.2

Dataset

GINO

CORAL

Transolver

GAOT

GAIA

Poisson-Gauss Elasticity 3D-Darcy

1.16 1.87 20.04

4.38 2.06 21.44

1.46 0.94 0.73

1.23 0.97 39.88

0.71 1.34 1.11

Boundary Value Problems

Figure 4 visualizes representative reconstructions on some MCB shapes. On the 3D Poisson BVP (Table 2), GAIA outperforms all baselines on every shape category, with the largest gains on gears (3.44% vs. 5.89% for NGF) and screws & bolts (3.68% vs. 5.52% for LNO). NGF is purpose-built for forward problems on linear symmetric PDEs; that GAIA exceeds it on every shape category, while remaining a general-purpose architecture, suggests the geometry-conditioned integral kernels carry meaningful advantage on real geometries. 4.2.3

Other Forward Problems

Table 3 reports forward-problem results on benchmarks where input and output share a discretization. GAIA achieves the lowest error on Poisson-Gauss (0.71%) and is competitive with the leading method on Elasticity (1.34% vs. 0.94% for Transolver). On 3D Darcy flow, Transolver achieves lower error than GAIA at the training resolution (0.73% vs. 1.11%). This advantage, however, is contingent on evaluation matching the training discretization. As shown in §4.3, both Transolver and GAOT degrade by an order of magnitude when evaluated at coarser resolutions on the Elasticity benchmark where they look strongest in-distribution; GAIA’s accuracy is stable. We view this as an architectural trade-off: GAIA’s integral-transform structure may prioritize discretization invariance over peak 8

in-distribution accuracy on shared-mesh forward tasks. For forward problems where deployment-time mesh resolution is fixed and matches training, Transolver remains a good choice; for applications where resolution varies—which includes many multi-query settings—GAIA’s robustness profile is preferable.

Figure 4: Results from the Mechanical Components Benchmark on the fitting, nut and gear categories. As shown in the error plots, GAIA reconstructs the solution with great fidelity across all shapes.

Median relative L 2 error (%)

Median Relative L 2 Error (%) Across Resolutions Median Relative L 2 Error (%)

60 50 40

Transolver GAIA GAOT

30 20 10 0

200

300

400

500

600

Resolution

700

800

900

8 6 4 2 0

972

Airfoil EIT RTE

10

0

1

2

3

4

5

6

7

8

9

10

Noise level η (% of ‖x‖∞)

(a) Discretization invariance on Elasticity (models trained at 972 points). GAIA maintains accuracy across (b) Graceful degradation in median relative L2 error resolutions while baselines degrade. (%) when increasing noise on inverse benchmarks.

Figure 5: Discretization invariance and noise robustness studies. 4.3

Discretization Invariance

A central claim of integral-transform operators like IAE-Net is discretization invariance (DI): accuracy must be maintained when the model is evaluated at point-cloud resolutions different from training. We test this on Elasticity—the forward benchmark where Transolver and GAOT achieve better results than GAIA—by training each model at the full resolution and evaluating on coarser subsets obtained via Farthest Point Sampling. Figure 5a shows that GAIA maintains accuracy across all tested resolutions, while Transolver and GAOT degrade significantly away from the training resolution. The in-distribution Transolver-vs-GAIA gap on Elasticity reported in Table 3 is an artifact of evaluating at the training resolution; under resolution mismatch, the comparison reverses. GAIA similarly outperforms Transolver at non-training resolutions on 3D Darcy (Appendix D.2, Figure 14). 4.4

Efficiency analysis

Table 4 reports inference cost on EIT (N = 2000 mesh nodes) on a single A6000 in FP32. Latency is the mean wall-clock time over 100 runs (batch size 1); throughput is measured at batch size 16. GAIA 9

Model

Params (M)

Latency (ms)

Memory (MB)

Throughput (samples/s)

GAIA 3.02 15.78 47.60 364.70 CORAL 0.84 4.99 32.70 252.9 GINO 2.31 36.54 1886.50 27.3 GAOT 3.40 10.45 45.4 135.9 Table 4: Computational cost comparison on EIT. GAIA has competitive latency and memory while maintaining the highest throughput among geometry-adaptive methods

achieves the highest throughput among geometry-adaptive methods (365 samples/s) at competitive latency (16 ms) and modest peak memory (48 MB). The closest method on throughput is CORAL at 253 samples/s. Among graph-based methods, GINO incurs roughly 2× GAIA’s latency at 40× the memory due to its dense FNO latent grid. 4.5

Noise Robustness

We evaluate GAIA’s robustness to test-time measurement noise by training once on clean data and evaluating at noise levels η ∈ {0%, 1%, 2%, 5%, 10%} using an additive Gaussian noise model scaled by the per-sample L∞ norm. Figure 5b shows graceful degradation across all three inverse benchmarks: at η = 10%, error grows by roughly 13× on Airfoil, 4× on EIT, and 7× on RTE relative to the clean baseline, with no sharp transitions or instability. The graceful-degradation behavior is consistent with what we would expect from an integral-transform architecture: noise in the input is averaged across the integration kernel rather than being amplified by point-wise operations. 4.6

Ablations

Table 5: Ablations. Top: Geometry conditioning: Cross-attention shows the best performance among conditioning methods. Bottom: Tokenizer pathways: Both tokenizers contribute to GAIA’s performance. Median relative L2 error (%). Conditioning Elasticity Airfoil Variant Elasticity Airfoil Concatenation FiLM Cross-attention

1.77 1.88 1.34

5.91 0.88 0.58

No boundary tokens No slice tokens Full model

2.28 1.90 1.34

1.06 1.03 0.58

Table 5 reports two ablations on Elasticity (forward) and Airfoil (inverse). Conditioning mechanism: cross-attention strictly outperforms both concatenation and FiLM, with the largest gap on Airfoil (0.58% vs. 5.91% for concatenation, an order of magnitude). This confirms that spatially adaptive conditioning—each query point attending independently to the geometry tokens—is essential for capturing local geometric features; a single global summary vector is insufficient on tasks with complex spatially varying physics. Tokenizer pathways: removing either the boundary tokens or the slice tokens degrades performance on both benchmarks. The two pathways encode complementary information: boundary tokens summarize global domain shape, while slice tokens capture spatially varying input physics. Sensitivity to token count and dimension is mild across the ranges we tested, with full sweeps in Tables (11-12) in Appendix E.

5

Conclusion and Limitations

GAIA conditions integral kernels on tokenized representations of domain geometry to solve forward and inverse problems on arbitrary domains in a single pass. It achieves state-of-the-art on all inverse and BVP benchmarks while maintaining stable accuracy across resolutions. Three limitations bound the present work. First, GAIA is designed for time-independent operators; extending the framework to time-dependent problems requires additional architectural treatment for temporal evolution and is left to future work. Second, on forward problems where input and output share a fixed dense discretization, GAIA trades a small amount of in-distribution accuracy for substantially improved 10

discretization invariance. Thirdly, the geometry-conditioning cross-attention scales as O(N K) in point-cloud size N and token count K; for very large 3D point clouds (N > 106 ), latent compression or sparse attention [50] may be needed to manage memory.

Datasets The new datasets used in our paper (beyond publicly available benchmarks) are available at https://drive.google.com/drive/folders/1OSrAGvJh14M5APg-oRvTKb0MQYQgS6VR? usp=sharing.

References [1] Roger Temam. Navier–Stokes equations: theory and numerical analysis, volume 343. American Mathematical Society, 2024. [2] Liliana Borcea. Electrical impedance tomography. Inverse problems, 18(6):R99–R136, 2002. [3] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed 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. [4] Sören Dittmer, Tobias Kluth, Peter Maass, and Daniel Otero Baguer. Regularization by architecture: A deep prior approach for inverse problems. Journal of Mathematical Imaging and Vision, 62(3):456–470, 2020. [5] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020. [6] Lu Lu, Pengzhan Jin, and George Em Karniadakis. DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193, 2019. [7] Shizheng Wen, Arsh Kumbhat, Levi Lingsch, Sepehr Mousavi, Yizhou Zhao, Praveen Chandrashekar, and Siddhartha Mishra. Geometry aware operator transformer as an efficient and accurate neural surrogate for PDEs on arbitrary domains. arXiv preprint arXiv:2505.18781, 2025. [8] Zongyi Li, Nikola Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, et al. Geometry-informed neural operator for large-scale 3D PDEs. Advances in Neural Information Processing Systems, 36:35836–35854, 2023. [9] Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: A fast transformer solver for PDEs on general geometries. arXiv preprint arXiv:2402.02366, 2024. [10] Yong Zheng Ong, Zuowei Shen, and Haizhao Yang. Integral autoencoder network for discretization-invariant learning. Journal of Machine Learning Research, 23(286):1–45, 2022. [11] Louis Serrano, Lise Le Boudec, Armand Kassaï Koupaï, Thomas X Wang, Yuan Yin, Jean-Noël Vittaut, and Patrick Gallinari. Operator learning with neural fields: Tackling PDEs on general geometries. Advances in Neural Information Processing Systems, 36:70581–70611, 2023. [12] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021. [13] Bing Yu et al. The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):1–12, 2018. [14] Justin Sirignano and Konstantinos Spiliopoulos. DGM: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339–1364, 2018. 11

[15] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced physicsinformed neural networks for forward and inverse PDE problems. Computer Methods in Applied Mechanics and Engineering, 393:114823, 2022. [16] Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, 43 (5):A3055–A3081, 2021. [17] Yinhao Zhu and Nicholas Zabaras. Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics, 366: 415–447, 2018. [18] Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Prediction of aerodynamic flow fields using convolutional neural networks. Computational Mechanics, 64(2):525–545, 2019. [19] Jonas Adler and Ozan Öktem. Solving ill-posed inverse problems using iterative deep neural networks. Inverse Problems, 33(12):124007, 2017. [20] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89):1–97, 2023. [21] Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized Fourier neural operators. arXiv preprint arXiv:2111.13802, 2021. [22] Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli. U-NO: U-shaped neural operators. arXiv preprint arXiv:2204.11127, 2022. [23] Jean Kossaifi, Nikola Kovachki, Kamyar Azizzadenesheli, and Anima Anandkumar. Multi-grid tensorized Fourier neural operator for high-resolution PDEs. arXiv preprint arXiv:2310.00120, 2023. [24] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485, 2020. [25] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755–6766, 2020. [26] Tapas Tripura and Souvik Chakraborty. Wavelet neural operator: a neural operator for parametric partial differential equations. arXiv preprint arXiv:2205.02191, 2022. [27] Tailin Wu, Takashi Maruyama, and Jure Leskovec. Learning to accelerate partial differential equations via latent global evolution. Advances in Neural Information Processing Systems, 35: 2240–2253, 2022. [28] Tailin Wu, Willie Neiswanger, Hongtao Zheng, Stefano Ermon, and Jure Leskovec. Uncertainty quantification for forward and inverse problems of PDEs via latent global evolution. In Proceedings of the AAAI conference on artificial intelligence, volume 38, pages 320–328, 2024. [29] Lianghao Cao, Thomas O’Leary-Roseberry, and Omar Ghattas. Derivative-informed neural operator acceleration of geometric MCMC for infinite-dimensional Bayesian inverse problems. Journal of Machine Learning Research, 26(78):1–68, 2025. [30] Sebastian Kaltenbach, Paris Perdikaris, and Phaedon-Stelios Koutsourelakis. Semi-supervised invertible neural operators for Bayesian inverse problems. Computational Mechanics, 72(3): 451–470, 2023. 12

[31] Aliaksandra Shysheya, Cristiana Diaconu, Federico Bergamin, Paris Perdikaris, José M Hernández-Lobato, Richard E Turner, and Emile Mathieu. On conditional diffusion models for PDE simulations. Advances in Neural Information Processing Systems, 37:23246–23300, 2024. [32] Roberto Molinaro, Yunan Yang, Björn Engquist, and Siddhartha Mishra. Neural inverse operators for solving PDE inverse problems. arXiv preprint arXiv:2301.11167, 2023. [33] Da Long, Zhitong Xu, Qiwei Yuan, Yin Yang, and Shandian Zhe. Invertible Fourier neural operators for tackling both forward and inverse problems. arXiv preprint arXiv:2402.11722, 2024. [34] Nicholas H Nelsen and Yunan Yang. Operator learning meets inverse problems: A probabilistic perspective. arXiv preprint arXiv:2508.20207, 2025. [35] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for PDEs on general geometries. Journal of Machine Learning Research, 24(388):1–26, 2023. [36] Sepehr Mousavi, Shizheng Wen, Levi Lingsch, Maximilian Herde, Bogdan Raonić, and Siddhartha Mishra. RIGNO: A graph-based framework for robust and accurate operator learning for PDEs on arbitrary domains. arXiv preprint arXiv:2501.19205, 2025. [37] Giovanni Catalani, Siddhant Agarwal, Xavier Bertrand, Frederic Tost, Michael Bauerheim, and Joseph Morlier. Neural fields for rapid aircraft aerodynamics simulations. Scientific Reports, 14 (1):25496, 2024. [38] Pan Du, Meet Hemant Parikh, Xiantao Fan, Xin-Yang Liu, and Jian-Xun Wang. Conditional neural field latent diffusion model for generating spatiotemporal turbulence. Nature Communications, 15(1):10416, 2024. [39] Tian Wang and Chuang Wang. Latent neural operator for solving forward and inverse PDE problems. Advances in Neural Information Processing Systems, 37:33085–33107, 2024. [40] Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying, Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. GNOT: A general neural operator transformer for operator learning. In International conference on machine learning, pages 12556–12569. PMLR, 2023. [41] Zijie Li, Kazem Meidani, and Amir Barati Farimani. Transformer for partial differential equations’ operator learning. arXiv preprint arXiv:2205.13671, 2022. [42] Zipeng Xiao, Zhongkai Hao, Bokai Lin, Zhijie Deng, and Hang Su. Improved operator learning by orthogonal attention. arXiv preprint arXiv:2310.12487, 2023. [43] Seungwoo Yoo, Kyeongmin Yeo, Jisung Hwang, and Minhyuk Sung. Neural Green’s functions. arXiv preprint arXiv:2511.01924, 2025. [44] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas. PointNet: Deep learning on point sets for 3D classification and segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017. [45] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. NeRF: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM, 65(1):99–106, 2021. [46] Alberto P Calderón. On an inverse boundary value problem. Computational & Applied Mathematics, 25:133–138, 2006. [47] Yuwei Fan and Lexing Ying. Solving optical tomography with deep learning. arXiv preprint arXiv:1910.04756, 2019. 13

[48] Sangpil Kim, Hyung-gun Chi, Xiao Hu, Qixing Huang, and Karthik Ramani. A large-scale annotated mechanical components benchmark for classification and retrieval tasks with deep neural networks. In European conference on computer vision, pages 175–191. Springer, 2020. [49] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012. [50] Angelos Katharopoulos, Apoorv Vyas, Nikolaos Pappas, and François Fleuret. Transformers are rnns: Fast autoregressive transformers with linear attention. In International conference on machine learning, pages 5156–5165. PMLR, 2020. [51] William RB Lionheart. EIT reconstruction algorithms: pitfalls, challenges and recent developments. Physiological measurement, 25(1):125, 2004. [52] Guillaume Bal, Alexandre Jollivet, and Vincent Jugnon. Inverse transport theory of photoacoustics. Inverse Problems, 26(2):025011, 2010. [53] Gerald Farin. Curves and surfaces for computer-aided geometric design: a practical guide. Elsevier, 2014. [54] Yixin Hu, Teseo Schneider, Bolun Wang, Denis Zorin, and Daniele Panozzo. Fast tetrahedral meshing in the wild. 2020. [55] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.

A

Datasets and Benchmarks

Below is a brief summary of the datasets that we use in the paper. A.1

Inverse Problem: Electrical Impedance Tomography

Electrical Impedance Tomography (EIT) is an imaging technique that involves solving an inverse problem to reconstruct the internal material distribution of an object using boundary measurements. It is governed by the electrostatic equations with variable conductivity given by the Laplace equation with mixed Dirichlet and Neumann boundary conditions: ∇ · (a(x)∇u(x)) = 0, for x ∈ Ω ⊆ Rd , d = 2, 3, .

(13)

Here, a(x) is the conductivity distribution of the medium. The goal is to recover this distribution using only the boundary data on ∂Ω using the Dirichlet-to-Neumann map (DtN map, also known as the Calderon operator [46]) Λa given by: Λa : H −1/2 (∂Ω) ∋ a∇u|∂Ω · n 7→ u|∂Ω ∈ H 1/2 (∂Ω). where H −1/2 (∂Ω) is the space of bounded linear functionals on the Sobolev space H 1/2 (∂Ω). In the case of full boundary measurement, it is a known result [46] that the inverse map F −1 : Λa 7→ a is well-posed and a is fully recoverable. However, in practice, we only have a finite amount of noisy data, making the EIT problem highly ill-posed and sensitive to noise [51]. We aim to directly learn the inverse map using boundary measurements across different geometries. The training and test datasets were generated by solving equation 13 on strictly star-shaped domains Ω defined by the polar radius r(θ) = r0 [1 + c1 cos(4θ) + c2 cos(8θ)]. Extending the data generation process in [32], the isotropic conductivity coefficient a(x, y) was drawn from a log-trigonometric probability distribution to ensure strict positivity and smoothness. Specifically, a(x, y) is constructed as: ! m X a(x, y) = exp ck sin(kπx) sin(kπy) , k=1

where m is drawn uniformly from {1, . . . , 5}, and the coefficients {ck }m k=1 are sampled independently from a uniform distribution U([−1, 1]). For each sampled conductivity, the governing equation 14

Table 6: Summary of benchmark problems. All datasets use median relative L2 error (%) as the evaluation metric. Benchmark

Governing Equations & Operator

Remarks

Poisson-Gauss [36]

∆u = f , x ∈ (0, 1)2 , u|∂Ω = 0. Operator: Gpg : L2 (Ω) → H01 (Ω) ∩ H 2 (Ω), f 7→ u. Source f is a superposition of M random Gaussians. ρ ∂tt u − ∇ · σ = 0. Operator: Gel : C(Ω) ∋ d 7→ σ ∈ L2 (Ω; R2×2 ), where d(x) = dist(x, Γinner ). Incompressible Rivlin–Saunders constitutive model. −∇ · (κ(x)∇u(x)) = 1, u|∂Ω = 0. Operator: Gdf : L∞ (Ω) → H01 (Ω), κ 7→ u. Lognormal permeability field; κ = exp(γ S/σS ).

Forward. Fixed unit-square domain. 128×128 uniform grid. 2056 train / 384 test samples.

Elasticity [8]

3D Darcy Flow

3D Poisson BVP [43]

−∆u = fanalytical , u|∂Ω = g(x, y, z). Oscillatory analytical source; randomized polynomial Dirichlet BCs.

EIT [32]

∇ · (a(x)∇u) = 0, mixed BCs. Operator: Λa : H −1/2 (∂Ω) → H 1/2 (∂Ω), a∇u · n 7→ u|∂Ω . Inverse map: Λa 7→ a(x).

Optical phy [47]

Airfoil [35]

Tomogra-

R v · ∇x u + σt u = σs Φ(v · v ′ )u dv ′ . Albedo 1 operator: Λ : L (∂Ω) → L1 (∂Ω), φ 7→ J+ . Inverse map: Albedo difference data ∆Λ 7→ δσs (x) scattering perturbation. Compressible Euler equations around a 2D airfoil (M∞ = 0.8, AoA = 0◦ ). Inverse map: sparse noisy Mach observations {ũ(xj )}j∈O 7→ u(x) over the full domain.

Forward. Unit square with randomly shaped interior hole (r ∈ [0.2, 0.4]). 972 unstructured nodes. 1,024 train / 256 test samples. Forward. Star-shaped 3D domains (spherical harmonic boundary, lmax = 4). Point cloud resolution 7000. 4096 train / 2048 test samples. BVP / Forward. 3D mechanical parts (MCB): gears, nuts, fittings, screws & bolts. Setting adapted from [43]. 200 shapes/category train, 20 test with 50 boundary conditions per shape. Inverse. Extends [32] to 2D starshaped domains (r(θ) = r0 [1 + c1 cos 4θ + c2 cos 8θ]). L = 20 Dirichlet patterns; M = 272 boundary sensors. Log-trigonometric conductivity. Inverse. Extends [47] to convex pentagon domains. S32 Discrete Ordinates solver; HG phase function (g = 0.9). 32 sources × 32 receivers. 6,500 samples. Inverse. Varying NACA airfoil geometries. 221 × 51 C-grid. 1% relative Gaussian noise, 10% sparse observations. 1,000 train / 100 test.

was solved numerically using quadratic Finite Element Method (FEM) under L = 20 distinct Dirichlet boundary conditions. The boundary data {gℓ }L ℓ=1 corresponds to plane waves incident from equispaced angles θℓ = 2π(ℓ−1) , formulated as: L u(x, y) ∂Ω = gℓ (x, y) = cos (2π(x cos θℓ + y sin θℓ )) . The corresponding Neumann measurement, Ψℓ = a ∂u ∂ν , representing the current flux, was evaluated at M = 272 boundary sensors for each excitation ℓ. Figure 6 shows the Dirichlet and Neumann measurements on the boundary for a given star-shaped domain and its corresponding recovered internal conductivity distribution. A.2

Inverse Problem: Optical Tomography

In optical tomography, the transport of photons through a scattering and absorbing medium occupying a bounded domain Ω ⊂ Rd is modeled by the stationary radiative transport equation (RTE). Let u(x, v) denote the particle density at spatial location x ∈ Ω and velocity direction v ∈ S d−1 . The medium is characterized by its scattering coefficient σs (x) ≥ 0 and absorption coefficient σa (x) ≥ 0. 15

(a) EIT: Dirichlet-Neumann boundary measurements.

(b) OT: Albedo matrix and scattering.

Figure 6: Inverse problems examples. (a) The EIT problem aims to reconstruct the internal conductivity distribution from boundary Dirichlet and Neumann values. (b) In OT, we reconstruct the scattering coefficient from the intensity at receivers given different source positions. The RTE is given by: Z

Φ(v · v ′ )u(x, v ′ )dv ′ ,

v · ∇x u(x, v) + σt (x)u(x, v) = σs (x)

x ∈ Ω,

(14)

S d−1

u(x, v) = ϕ(x, v),

(x, v) ∈ Γ− .

(15) ′

Here, σt (x) = σa (x) + σs (x) is R the total attenuation coefficient, and Φ(v · v ) is the scattering phase function, normalized such that S d−1 Φ(v · v ′ )dv ′ = 1. The boundary phase space is partitioned into inflow (Γ− ) and outflow (Γ+ ) boundaries, defined as Γ± = {(x, v) ∈ ∂Ω × S d−1 : ±nx · v > 0},

(16)

where nx is the unit outer normal vector at x ∈ ∂Ω. In our experiment, following [47], both the boundary illumination and the recorded data are assumed to be independent of the velocity variable. Specifically, the prescribed incoming particle density is modeled as an isotropic source depending only on the spatial variable, ϕ(x) ∈ L1 (∂Ω). Similarly, the measurements on the outflow boundary record the outgoing angular flux (or current), integrating the particle density over all outward-pointing R velocities: J+ (x) = v·nx >0 (v · nx )u(x, v)dv. The associated boundary observation operator, known as the Albedo operator Λ, thus maps the spatial input on the boundary to the spatial output: Λ : L1 (∂Ω) → L1 (∂Ω),

Λ : ϕ(x) 7→ J+ (x).

(17)

The continuous inverse problem seeks to recover the unknown medium properties, characterized by the scattering and absorption coefficients, from the knowledge of Λ. Formally, this defines the continuous inverse map F −1 from the space of bounded linear operators to the space of continuous medium parameters:
F −1 : L L1 (∂Ω), L1 (∂Ω) → C(Ω) × C(Ω), F −1 (Λ) = (σs , σa ). The well-posedness and Lipschitz stability of this inverse map have been established under suitable conditions [52]. For the numerical experiments and data generation in this work, we restrict our domain Ω ⊂ R2 to randomly scaled convex pentagons. We assume the absorption coefficient is a known, spatially invariant constant, σa (x) ≡ 0.01, reducing the inverse problem to the recovery of the spatially varying scattering coefficient σs (x). The highly forward-peaked anisotropic scattering in biological tissues is modeled using the 2D Henyey-Greenstein (HG) phase function, Φ(v · v ′ ) =

1 − g2 , 2π(1 + g 2 − 2g(v · v ′ ))

where the anisotropy factor is set to g = 0.9. The spatial distribution of the scattering coefficient is constructed using a randomized superposition of low-frequency trigonometric modes passed through a scaled logistic sigmoid function. Let x′ = (x′1 , x′2 ) denote the spatial coordinates normalized by the domain scale factor. The scattering field is defined as: !!−1 m X σs (x) = 1.0 + 1 + exp −3 ck sin(ωx,k πx′1 ) sin(ωy,k πx′2 ) , k=1

16

where the number of modes m is drawn uniformly from {2, . . . , 5}, the expansion coefficients ck are sampled independently from a uniform distribution U([−1, 1]), and the spatial frequencies ωx,k and ωy,k are drawn uniformly from the integer set {1, 2, 3}. This construction ensures that the scattering coefficients are smooth and strictly bounded within the physical range σs (x) ∈ (1.0, 2.0). The continuous boundary illumination is discretized using Nsrc = 32 localized isotropic point sources distributed uniformly along ∂Ω. The forward transport problem is solved using the Discrete Ordinates (S32 ) method with angular discretization of 32 angles, and a Finite Element spatial discretization, using GMRES for the source iteration. The corresponding outgoing current is measured at Nrec = 32 interleaved receiver locations. This discrete setup approximates the Albedo operator as a dense matrix Λ ∈ R32×32 . To isolate the scattering anomalies and suppress structural geometric artifacts, we employ a difference imaging approach. A baseline Albedo matrix Λ0 is simulated using a known, constant background scattering coefficient σs,0 = 1.0. The objective of the discrete inverse problem is then to infer the scattering perturbation δσs (x) = σs (x) − σs,0 strictly from the measured difference data ∆Λ = Λ − Λ0 .

Figure 7: Model predictions on RTE. Comparison of ground truth solution against predictions from GAIA, GAOT, Transolver, CORAL, and GINO.

A.3

Inverse Problem: Airfoil Flow Reconstruction

We consider the problem of reconstructing a full aerodynamic flow field from sparse, noisy measurements over a family of transonic solutions over airfoil geometries. The underlying physical system is governed by the compressible Euler equations, ∂ρf + ∇ · (ρf v) = 0, ∂t

∂ρf v + ∇ · (ρf v ⊗ v + pI) = 0, ∂t

∂E + ∇ · (E + p)v = 0, (18) ∂t

where ρf is the fluid density, v is the velocity vector, p is the pressure, and E is the total energy. Viscous effects are neglected. Far-field conditions are prescribed as ρ∞ = 1, p∞ = 1.0, M∞ = 0.8, AoA = 0◦ , with a no-penetration condition imposed at the airfoil surface. Here M∞ is the Mach number and AoA is the angle of attack. The dataset is drawn from the NACA airfoil corpus of [35], comprising N = 1,200 flow solutions, one per airfoil shape, generated with a second-order implicit finite volume solver on a body-fitted C-grid of Nx × Ny = 221 × 51 = 11,271 quadrilateral elements, refined near the airfoil surface. Airfoil geometries are parameterized via a design element approach [53] in which the baseline NACA0012 profile is deformed by displacing the control nodes of an enclosing cubic design element in the vertical direction, with displacements drawn from d ∼ U [−0.05, 0.05]. The dataset is partitioned into Ntr = 1,000 training, Nval = 100 validation, and Nte = 100 test configurations. We formulate it as a sparse-to-full field reconstruction problem. The quantity of interest is the steadystate Mach number field u : Ω → R, discretized over all 11,271 mesh nodes. At inference time, only a sparse subset of observations {u(xj )}j∈O , with O ⊂ {1, . . . , Nx Ny }, is available as input, and these observations are further corrupted by zero-mean Gaussian noise scaled to the per-sample peak Mach number:
ũi (xj ) = ui (xj ) + ηi,j , ηi,j ∼ N 0, σ 2 ∥ui ∥2∞ , j ∈ O, (19) with σ = 0.01, corresponding to 1% of the per-sample L∞ norm. The task is to recover the complete Mach field u at all grid nodes, including those with no sensor coverage. The model must simultaneously denoise, interpolate, and extrapolate from incomplete observations over geometrically varying domains. 17

Figure 8: Model predictions on inverse benchmarks. Top: EIT conductivity reconstruction. Bottom: Airfoil Mach field reconstruction from sparse noisy observations. Ground truth is shown on the left, followed by predictions from baselines. GAIA solutions show good agreement with the ground truth solutions while other models struggle. A.4

3D Boundary Value Problem: Poisson on MCB

To test the generalization of our framework to unseen problem domains, we construct a new PDE dataset using the Mechanical Components Benchmark (MCB) dataset [48], which contains a variety of 3D mechanical part shapes. We generate tetrahedral meshes for shapes in four categories (SCREWS & BOLTS, NUTS, FITTINGS, and GEARS) by meshing the interiors of unit-cube-normalized shapes using fTetWild [54]. The shape collection is divided into 200 shapes for training and 20 shapes for testing per category. Shapes within the same category exhibit diverse geometries, presenting a standard challenge for learned solution operators to generalize effectively. We define a highly non-linear, oscillatory analytical source term fanalytical (x, y, z). The base function is constructed using trigonometric variations along all three spatial axes: fanalytical = T1 (x, y) + T2 (x, y) + T3 (z) where the individual components are defined by the constants A = 1.25, B = 1.5, C = 1.5, and D = 1.5, such that:
T1 (x, y) = − (Aπ)2 + (ACπ)2 sin(Aπx) cos(ACπy) T2 (x, y) = (Aπ)2 cos(Aπx) (1 − sin(ABπy)) + (ABπ)2 sin(ABπy) (1 − cos(Aπx)) T3 (z) = 2(ADπ)2 cos(2ADπz) To generate multiple unique samples per geometry, we randomize the Dirichlet boundary conditions g(x, y, z) at the mesh surface. The boundary values are prescribed using a parameterized polynomial function: g(x, y, z) = E(x3 − 3xy 2 ) + F (y 3 − 3x2 y) + (x2 − z 2 ) For each generated sample, the coefficients E and F are sampled from uniform distributions such that E ∼ U(−1.0, 1.0) and F ∼ U (0.0, 1.0). To generate the dataset, we employ the FEniCS finite element framework [49] to solve Poisson’s equation on the generated meshes. A.5

Forward Problem: 3D Darcy Flow

We consider the stationary Darcy flow through a heterogeneous porous medium. Let the spatial domain Ω ⊂ R3 represent three-dimensional domain. The fluid pressure u(x) is modeled by the elliptic partial differential equation: −∇ · (κ(x)∇u(x)) = 1,

x ∈ Ω,

subject to homogeneous Dirichlet boundary conditions, u|∂Ω = 0. Here, κ(x) represents the strictly positive permeability field of the medium. The operator learning objective is to approximate the 18

Figure 9: Boundary Values on MCB shapes

Figure 10: Results on the screws&bolts category from the Mechanical Components Benchmark.

continuous solution operator Gd that maps the permeability field, defined on a variable geometries, to the resulting scalar pressure field. Formally, for a given domain Ω, the mapping is defined as: Gdf : L∞ (Ω) → H01 (Ω),

Gdf (κ) = u.

To learn this operator, we construct a dataset by solving the Darcy flow equation on Ω ⊂ R3 . The boundary ∂Ω of each sample is generated as a "star shape" defined by a random spherical harmonic expansion up to degree lmax = 4: r(θ, ϕ) =

4 l X X

clm Ylm (θ, ϕ),

clm ∼ N (0, σl ).

(20)

l=0 m=−l

The resulting geometries are discretized into watertight unstructured tetrahedral meshes. The permeability field κ(x) is computed as a continuous log-normal random field, κ(x) = exp(γ · S(x)/σS ), where S(x) is a superposition of M = 30 random Fourier modes: S(x) =

M X

αj cos(kj · x + ϕj ).

(21)

j=1

Here, the wave vectors kj ∼ N (0, σk2 I) control the spatial correlation length, phases ϕj ∼ U[0, 2π] ensure translation invariance, and amplitudes αj ∼ N (0, 1) provide variance. The field is normalized by its sample standard deviation σS and scaled by γ = 0.5 to ensure a smooth, strictly positive, and spatially heterogeneous medium. The forward problem is solved using the standard Galerkin Finite Element Method (FEM) implemented in FEniCS. We utilize linear Lagrange (P1 ) basis functions for both the solution space and the coefficient field. Following the FEM solution, we extract the training data by randomly sampling Np = 7000 points within the domain volume for each simulation. A.6

Forward Problem: Elasticity Problem

We evaluate the model on a standard Elasticity benchmark, where the objective is to learn the mapping from a domain’s geometry to its resulting internal stress field under external loading. The general force balance of a solid body is governed by the hyperelastic equation: ρ∂tt u − ∇ · σ = 0, 19

(a) Elasticity

(b) Poisson-Gauss

Figure 11: Input output fields where ρ represents the mass density, u the displacement field, and σ the Cauchy stress tensor. The strain field is related to the displacement field via standard kinematic relations, and the system is closed by a hyperelastic constitutive model describing the non-linear stress-strain relationship. For this benchmark from [35], the stationary solution is considered for a unit square hyper-elastic incompressible Rivlin-Saunders specimen with a hole at its center. The geometry of the hole is randomly sampled such that the radius always takes a value between 0.2 and 0.4 (see [35] for details). The specimen is subjected to mixed boundary conditions: it is clamped at the bottom boundary (u = 0) and is under a constant vertical tension traction on its top boundary. The stress field is computed for different geometries with a finite elements solver with about 100 quadratic quadrilateral elements. The target function in this dataset is the stress field σ, which is available at 972 unstructured coordinates. The operator learning task thus seeks to approximate the continuous solution operator Gel that maps the geometry of the domain directly to the stress field. The geometry is encoded as the shortest-path distance function from the inner boundary of the hole, d(x) = dist(x, Γinner ). The formal operator mapping is thus defined as: Gel : C(Ω) ∋ d → σ ∈ L2 (Ω; R2×2 ). Figure 11a demonstrates an example of a sample with a deformation in the center. The input field is distance of the unordered point cloud from the inner boundary.

Figure 12: Model predictions on Elasticity. Comparison of ground truth solution against predictions from GAIA, GAOT, Transolver, CORAL, and GINO. A.7

Forward Problem: Poisson-Gauss

To demonstrate GAIA’s performance on a standard fixed-grid benchmark, we consider the classical Poisson’s equation [36]. Let the domain be the unit square Ω = (0, 1)2 . The governing PDE is given by: ∆u(x) = f (x), x ∈ (0, 1)2 , (22) subject to homogeneous Dirichlet boundary conditions, u|∂Ω = 0. The objective of the neural operator is to learn the inverse Laplacian operator Gpg = ∆−1 , which maps the right-hand side source function to the scalar solution field. Formally, this defines the continuous solution operator: Gpg : L2 (Ω) → H01 (Ω) ∩ H 2 (Ω),

Gpg (f ) = u.

(23)

The source term f is procedurally generated as a superposition of M random Gaussian pulses:   M X (x − µx,i )2 + (y − µy,i )2 f (x, y) = exp − , (24) 2σi2 i=1 20

Figure 13: Model predictions on Poisson-Gauss. Comparison of ground truth solution against predictions from GAIA, GAOT, Transolver, CORAL, and GINO. where M is an integer drawn from a geometric distribution. The mean coordinates are sampled uniformly such that µx,i , µy,i ∼ U(0, 1), and the standard deviations are sampled as σi ∼ U(0.025, 0.1). The ground truth data pairs (f, u) are evaluated on a high-resolution 128 × 128 uniform grid. We refer to [36] for further dataset details. Figure 11b visualizes a representative input source function and its corresponding computed solution.

B

GAIA Implementation Details

All models were trained on a single NVIDIA A6000 GPU using the Adam optimizer [55] with cosine decay learning rate scheduling and 5% warmup starting at 10% of the peak learning rate. We use 8 attention heads, Fourier positional encodings with 10 frequency bands and base frequency ω0 = 30. All hyperparameters for experiments are given in Table 7. Table 7: Training and architectural hyperparameters across benchmarks. BS refers to batch size, ep.: epochs, Nin /Nout : number of input/output points per sample, GHD: Hidden dimension, GAD: Token dimension x Number of tokens, AD: attention dimension, NGT: Number of geometric tokens. Dataset

Ntrain Ntest

OT EIT Airfoil Elasticity P-Gauss 3D Darcy

8000 4096 1000 1024 2048 4096

MCB: Gear MCB: Screw MCB: Fitting MCB: Nut

C

10k

Nin

Training Nout

2000 32 × 32 ∼5000 2048 272 2000 100 1127 11271 256 972 972 256 1282 1282 2048 7000 7000 1k

∼6700 ∼15800

Architecture Geometry Conditioning BS LR Ep. Modes Width Blocks GHD GAD AD NGT. 4 4

2e-4 1e-4

8

5e-4

500

96 32 96 128

2e-4 1

2e-4 1e-4 100 1e-4 2e-4

128 96

128

4 5 4

64 128

5

128

4

256

256 128

8 8 4 16 8 8

128

128

4

256

256

128

64 32 128 256 128

64

256

Baselines

Unless specifically noted, all baseline implementations follow their default hyperparameters in their published repositories. While CORAL and LNO can technically handle varying input/output spatial domains, the codebase had to be modified slightly to accomplish this. C.1

CORAL

CORAL [11] casts operator learning as a regression in the latent spaces of two implicit neural representations (INRs). Each input sample ai and output sample ui is encoded by fitting a shiftmodulated SIREN via auto-decoding: a shared base network is trained jointly with per-sample latent codes zai , zui , where the codes are passed through small hypernetworks that produce additive shifts on the SIREN activations. The codes are obtained through a short inner loop of gradient descent on the reconstruction loss. Once the two INRs are trained, a separate network learns the mapping zai 7→ zui in latent space. Because the INRs are queried point-wise, CORAL natively supports varying input and output meshes across samples. 21

8

We use the official implementation3 . Both INRs are 4-layer SIRENs of hidden width 256 with shift-only modulation generated by a single-layer hypernetwork of width 128. The latent regression network za 7→ zu is a 3-block residual MLP with Swish activations, trained with Adam (no weight decay) and a step decay of γ = 0.9. The two-stage protocol trains the INRs first, then freezes them and fits the regression network. We meta-learn the inner learning rate at a meta-rate of 10−4 . As the official codebase assumes a shared mesh between input and output, we modified the data pipeline to query the two INRs on distinct point clouds for boundary-value problems and inverse problems where sensor and reconstruction grids differ. C.2

GINO

GINO [8] is a hybrid neural operator that combines a graph neural operator (GNO) [24] with a Fourier neural operator (FNO) to handle inputs and outputs on irregular point clouds. An input GNO encoder integrates the input function over local neighborhoods of each latent grid point, lifting the irregular sample onto a regular Cartesian latent grid. An FNO processor operates on this latent grid, exploiting the FFT for global integration. An output GNO decoder then projects the processed latent field back onto an arbitrary query point cloud. The encoder and decoder rely on radius graphs: each output node attends to all input nodes within a fixed radius r. We use the official implementation4 . The FNO processor uses 6 layers with hidden width 32, group normalization, and a channel MLP expansion of 0.5; AdaIN conditioning is enabled with 8 feature channels. Both GNO encoder and decoder use linear kernel transforms with a half-cosine kernel weighting. The input GNO MLP has hidden layers [128, 256, 128] and the output GNO MLP has hidden layers [256, 512, 256]. The GNO radius r is tuned to the characteristic length scale of each domain and is shared between the input and output GNO; values are listed in Table 8. Table 8: GINO GNO radius r across benchmarks. The same radius is used for the input and output GNO.

C.3

Dataset

EIT

Airfoil

Poisson-Gauss

Elasticity

RTE

3D Darcy

r

0.25

0.125

0.125

0.125

5.0

0.25

Transolver

Transolver [9] is a transformer-based neural operator that adapts standard self-attention to PDE solving on irregular meshes. Rather than computing attention directly between mesh nodes – which scales quadratically with the point-cloud size – it introduces a Physics-Attention mechanism that (softly) assigns each mesh point to a small set of learnable physical states (or slices), runs self-attention over the slice representations, and then projects back to the original mesh. This reduces attention cost from O(N 2 ) to O(N K) for K slices and lets the model scale to large unstructured meshes. Because both the slice assignment and the back-projection share the same point set, Transolver requires the input and output to live on a common discretization, which is why we omit it from inverse and BVP comparisons. We use the official implementation5 . The model uses 8 Transolver blocks of hidden width 256 with 8 attention heads, an MLP expansion ratio of 2, and 32 slice tokens per block; the slice projection uses 8 reference points. Training uses AdamW with a learning rate of 10−3 , weight decay 10−5 , and batch size 8. C.4

Geometry Aware Operator Transformer (GAOT)

GAOT [7] is a graph-based transformer architecture for operator learning on unstructured grids. It uses a multiscale graph encoder (MAGNO) to aggregate point-cloud features onto a fixed Cartesian latent token grid, applies transformer attention over the latent tokens, and decodes back to the mesh through a corresponding graph decoder. For the forward problems (Elasticity, Poisson-Gauss and 3D Darcy), we adopt the default hyperparameters provided in the GAOT repository without modification 6 . For 2D problems, the latent token 3 5

4 https://github.com/neuraloperator/neuraloperator https://github.com/LouisSerrano/coral 6 https://github.com/thuml/Transolver https://github.com/camlab-ethz/GAOT

22

grid is 64×64; for the 3D Poisson problem, we use a 64×32×32 latent grid as specified in the original paper [7]. Training uses the default AdamW optimizer with a mixed learning rate schedule (linear warmup followed by cosine decay). For the three inverse problems (EIT, airfoil flow and OT), we use a modified version of the GAOT codebase. Since inverse problems involve distinct input and output domains (e.g., boundary measurements to interior fields), we modify the MAGNO graph construction to build separate encoder and decoder neighborhood graphs, rather than a single shared graph as in the original implementation. We additionally lower the initial learning rate to 1e−4 , as the default 8e−4 did not yield good performance on this problem. All other architectural hyperparameters remain at their defaults. C.5

Latent Neural Operator (LNO)

LNO [39] learns operator mappings via attention through a low-dimensional latent space. The encoder applies cross-attention between a fixed set of learnable latent queries and the input point cloud, producing a compact latent representation. A stack of self-attention blocks then processes these latent tokens. Finally, the decoder applies cross-attention from the output query points to the processed latent tokens to produce the predicted field. We apply LNO to the four forward Poisson problems on MCB dataset, using the Plasticity config from the original repository7 as the reference. The architectural hyperparameters are kept identical. Batch size is reduced to 1 to accommodate the MCB meshes in memory, and the model is trained for 100 epochs like our models. For the inverse problems, we use the LNO propagator config from the original repository as the reference, with batch size set to 64 for speed. C.6

Neural Greens Function (NGF)

NGF [43] approximates the solution operator of linear symmetric PDEs by learning a low-rank decomposition of the Green’s function G(x, y) on a point cloud. Two coordinate networks predict left- and right-vector bases at each point, and the solution at a query is recovered by integrating the learned Green’s function against the source field. The architecture is designed for forward problems whose governing operator is linear, self-adjoint, and admits a Green’s function representation. It is not applicable to inverse problems or nonlinear forward problems. So, we apply NGF to four Poisson problems on irregular geometries (fitting domain, gear, nut, screws/bolts) using the default network and training configuration from the original repository8 without modification.

D

Additional Results

D.1

Noise Robustness

We adopt a relative-Gaussian noise model scaled by the per-sample L∞ norm of the measurement field. Let x denote a clean measurement vector; the corrupted measurement is: 
2  x̃ = x + ε, ε ∼ N 0, η ∥x∥∞ I . (25) for noise levels η ∈ {0%, 1%, 2%, 5%, 10%} The three benchmarks differ in what the measurement vector x is and at what granularity the L∞ norm is computed. • Airfoil. The measurement is the partial Mach field observed at a fraction of mesh points. A single ∥x∥∞ is computed per sample over all observed points. • EIT. The measurement consists of L = 20 paired Dirichlet–Neumann boundary patterns evaluated at M = 272 sensors. • RTE. The measurement is the difference Albedo matrix ∆Λ ∈ RNrec ×Nsrc . Because receivers have very different outgoing flux magnitudes, we use a per-receiver-row L∞ scale: εr,s ∼ N (0, (η ∥∆Λr,: ∥∞ )2 ). 7

https://github.com/L-I-M-I-T/LatentNeuralOperator

23

8

https://github.com/KAIST-Visual-AI-Group/NGF.git

For each benchmark and noise level, we evaluate the trained model under 10 independent noise realizations and report the mean and standard deviation of the median relative L2 error across these seeds; results are summarized in Table 9. Across all three benchmarks, GAIA degrades gracefully under measurement noise. Seed-to-seed variation remains under 3% of the mean at every noise level, showing that this robustness is a stable property of the trained model. Table 9: Median relative L2 error (%) under increasing noise levels η on the Airfoil, EIT, and RTE inverse benchmarks. Values are reported as mean ± std over 10 independent noise realizations. Noise level η Benchmark Airfoil EIT RTE

D.2

0%

1%

2%

5%

10%

0.49 ± 0.03 0.70 ± 0.00 1.42 ± 0.00

0.59 ± 0.03 0.76 ± 0.01 1.53 ± 0.01

1.12 ± 0.06 0.91 ± 0.01 1.83 ± 0.01

3.87 ± 0.14 1.53 ± 0.03 3.57 ± 0.03

6.51 ± 0.09 2.79 ± 0.04 10.16 ± 0.09

Discretization invariance for Transolver

We conduct a similar discretization invariance study for 3D Darcy benchmark (see Figure 14). While Transolver’s degradation is less severe than for the Elasticity benchmark, GAIA still outperforms Transolver at non-training resolutions.

Figure 14: Discretization invariance on 3D Darcy. Although Transolver achieves lower error at the training resolution (100%), GAIA degrades more gradually under subsampling and achieves lower error at coarser resolutions.

D.3

Data augmentation for discretization invariance

Discretization invariance can be further improved by training on multiple resolutions simultaneously. We demonstrate this on the Poisson-Gauss benchmark, comparing models trained with and without multi-resolution data augmentation. The data augmentation used in IAE-Net is implemented here via a stride-based multi-resolution training scheme. The original data is generated on a 128 × 128 grid and subsampled to 64 × 64 for standard training. For augmented training, the 64 × 64 data is upsampled back to the canonical 128 × 128 resolution using bicubic interpolation. Multi-resolution views are then generated by striding this grid. Each training iteration simultaneously optimizes over multiple resolutions R = {32 × 32, 32 × 64, 64 × 32, 64 × 64, 64 × 128, 128 × 64, 128 × 128}, with the total loss: X Ltotal = λT · L(ûT , uT ), (26) T ∈R

24

0.30

With DA Without DA

Mean Relative L 2 Error

0.20 0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03

8 12 8x

64 8x

28 x1

12

12

Test Resolution

64

4 x6 64

2 x3 64

4 x6 32

32

x3

2

0.02

Figure 15: Discretization invariance on Poisson-Gauss. Multi-resolution augmentation (blue) provides one-shot generalization across resolutions compared to single-resolution training (orange).

where λT are per-resolution weights (set to 1.0 in this experiment). As shown in Figure 15, multiresolution training yields stable performance across all tested resolutions. We note that for unstructured point clouds, generating equivalent multi-resolution training views would require scattered data interpolation to obtain field values at new point locations, which is less straightforward than grid striding and may introduce interpolation artifacts. D.4

Efficiency analysis

Table 4 reports computational costs on the EIT benchmark (N = 2000 mesh nodes) measured on a single NVIDIA A6000 GPU in FP32. Latency is the mean wall-clock time over 100 runs (batch size 1) after 20 warm-up iterations with CUDA synchronization; throughput is measured at batch size 16. GAIA achieves competitive latency and memory while maintaining the highest throughput among geometry-adaptive methods. To assess the computational efficiency of GAIA relative to baseline methods, we conduct a comprehensive timing analysis on the EIT benchmark with N = 2000 mesh nodes. All experiments are performed on a single NVIDIA A6000 GPU using FP32 precision. Metrics. We evaluate the following computational metrics: • Parameters: Total number of learnable weights in the model, reported in millions (M). • Latency: Wall-clock time (ms) for a single inference with batch size 1. We report the mean over 100 timed runs following 20 warm-up iterations to ensure stable GPU clock speeds and memory allocation. CUDA synchronization is enforced before each timing measurement to account for asynchronous kernel execution. • Peak Memory: Maximum GPU memory allocated during a forward pass with batch size 1, measured after resetting CUDA memory statistics. • Throughput: Maximum number of samples processed per second with batch size 16, measured by timing batched inference over 50 iterations. D.5

Training Randomness

To assess sensitivity to training randomness (weight initialization, data shuffling), we train GAIA five times with different random seeds on the EIT and Airfoil benchmarks. Table 10 reports the resulting statistics. The small standard deviation indicates stable performance across runs. For airfoil, we remove the random noise to isolate model randomness. 25

Table 10: Seed stability on EIT and airfoil datasets. Median relative L2 error (%) over 5 independent training runs. Dataset Error [%] (Mean ± Std) 0.71 ± 0.04 0.49 ± 0.03

EIT Airfoil

E

Ablations

We examine the sensitivity of the model to the tokenization hyperparameters, and provide details on the choice of conditioning mechanism, and the contribution of each token pathway. First, we conduct ablation studies on the EIT benchmark to investigate the sensitivity of GAIA to the token dimension and the total number of tokens of the tokenization module. We additionally ablate the geometry conditioning mechanism itself. E.1

Token Dimension

Table 11 reports the median relative L2 error as the token dimension varies from 4 to 64, with the number of tokens fixed at 8. Performance is best at 4 and degrades gradually as the dimension increases, suggesting that a compact token representation is sufficient to capture the relevant geometric information for this benchmark. Larger token dimensions may introduce unnecessary capacity that hinders optimization without providing additional representational benefit. Token Dim

4

8

16

32

64

L2 Median 0.00698 0.00712 0.00737 0.00804 0.00832 Table 11: Effect of token dimension with number of tokens fixed to 8.

E.2

Number of Geometric Tokens

Table 12 examines the effect of varying the token count from 4 to 64 with the token dimension fixed at 8. Performance improves sharply from 4 to 8 tokens, after which it plateaus. This indicates that for simple shapes, modest number of tokens is sufficient to encode the geometric structure of the domain, and that additional tokens do not meaningfully improve accuracy. Num Tokens

4

8

16

64

L2 Median 0.00851 0.00712 0.00727 0.00760 Table 12: Effect of number of tokens with token dimension fixed to 8.

E.3

Number of Slice Tokens

The slice tokenizer aggregates the input sensor data into a fixed set of physics-informed tokens via softmax-weighted clustering. Table 13 examines the effect of varying the slice-token count from 32 to 96, with all other hyperparameters held fixed. Performance is best at 64 slice tokens and degrades only mildly outside this setting, indicating that a moderate number of slices is sufficient to summarize the boundary measurements without introducing redundancy. E.4

Number of Blocks

Table 14 examines the effect of varying the number of blocks from 3 to 6 on the EIT benchmark, with all other hyperparameters held fixed. Performance is best at 5 blocks and remains close at neighbouring depths, indicating that the model is robust to this choice. 26

Num Slice Tokens

32

48

64

96

L2 Median 0.00797 0.00873 0.00712 0.00735 Table 13: Effect of the number of slice tokens on the EIT benchmark, with all other hyperparameters fixed. Num Blocks

3

4

5

6

7

L2 Median 0.00762 0.00789 0.00712 0.00739 0.008406 Table 14: Effect of number of blocks with all other hyperparameters fixed.

E.5

Geometry Conditioning Ablation Details

Table 5 compares three alternatives for using geometry information to condition the integral kernel on a forward problem (Elasticity) and an inverse problem (Airfoil). In concatenation, the geometry tokens are mean-pooled into a single vector, projected, and concatenated onto the per-point feature before entering the kernel MLP, thus providing global context but no spatial selectivity. In FiLM conditioning, the mean-pooled token vector is projected to produce per-channel scale and shift parameters that affinely modulate the combined feature representation. The proposed cross-attention mechanism instead lets each spatial query point attend independently to the full set of geometry tokens via multi-head attention. Cross-attention consistently achieves the lowest error on both benchmarks, implying that spatially adaptive conditioning is important for capturing local geometric features. E.6

Tokenizer Ablation Details

The IAE refinement stage stacks several adaptive blocks, each applying integral transforms followed by token-updated cross-attention. Table 5 evaluates the contribution of each token pathway by removing either the boundary feature tokens (produced by the PointNet-style encoder from Xbnd ) or the slice tokens (produced by the soft-clustering encoder from the interior mesh and field values). Both ablations increase error relative to the full model on the Elasticity and Airfoil benchmarks, implying that both tokenization pathways contribute information: boundary tokens capture the global domain shape, while slice tokens encode the spatially varying physics of the input field.

27