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 [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 [7], [8], 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 geometry-adaptive integral operator that conditions learned kernels on domain geometry via cross-attention, 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 \(L^2\) 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 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.1 Problem Formulation↩︎

We consider a generic time-independent PDE, \[\mathcal{D}(c, u) = q, \quad \forall x \in D \subset \mathbb{R}^d, \quad \mathcal{B}(u) = u_b, \quad x \in \partial D,\] where \(u : D \to \mathbb{R}^m\) is the solution, \(c \in \mathcal{C}\) represents the physical parameters, \(q \in \mathcal{F}\) is the forcing term, \(u_b\) prescribes the boundary values; and \(\mathcal{D}\) and \(\mathcal{B}\) are the differential and boundary operators, respectively. Also, let \(\chi_D\) denote the domain indicator. Here \(\mathcal{U}\), \(\mathcal{C}\), and \(\mathcal{F}\) denote suitable function spaces on \(D\). Encapsulating the problem setup into \(a = (c, q, u_b, \chi_D) \in \mathcal{A}\), the forward operator \(\mathcal{S} : \mathcal{A} \to \mathcal{U}\) maps \(a \mapsto 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 \(\mathcal{M} : \mathcal{U} \to \mathcal{V}\) mapping \(u\) to an observation \(v = \mathcal{M}(u)\). The generalized forward operator \(\mathcal{F} = \mathcal{M} \circ \mathcal{S} : \mathcal{A} \to \mathcal{V}\) maps parameters directly to observables, and the inverse operator \(\mathcal{F}^\dagger : \mathcal{V} \to \mathcal{C}\) recovers \(c = \mathcal{F}^{\dagger}(v, \chi_D)\) from these measurements.

Throughout, let \(\Psi : \mathcal{X} \to \mathcal{Y}\) denote the target mapping (\(\Psi \equiv \mathcal{S}\) in the forward case, \(\Psi \equiv \mathcal{F}^\dagger\) in the inverse case). Our goal is to approximate \(\Psi\) given access to data pairs \((f^{(i)}, g^{(i)})\) with \(g^{(i)} = \Psi(f^{(i)})\).

3.2 IAE-Net↩︎

IAE-Net [10] learns an operator \(\Psi : \mathcal{X} \to \mathcal{Y}\) between function spaces on compact domains \(\Omega_x \subset \mathbb{R}^{d_x}\) and \(\Omega_y \subset \mathbb{R}^{d_y}\), from finite discretizations \(S_x \subset \Omega_x\), \(S_y \subset \Omega_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 \to a_0 \to a_1 \to \dots \to a_L \to g\), with each intermediate function \(a_i\) defined on \(\Omega_a = [0,1]^{d_a}\). Each block performs an encode–process–decode transformation \(a \to v \to u \to b\), where \(v\) and \(u\) live on a fixed latent domain \(\Omega_z = [0,1]^d\), discretized on a grid \(S_z = \{z_j\}_{j=1}^m\) independent of the input: \[\small v(z) = \int_{\Omega_x} \phi_1(a(x), x, z; \theta_1)\, a(x)\, dx,\;u(z) = \phi_0(v(z); \theta_0),\;b(y) = \int_{\Omega_z} \phi_2(u(z), y, z; \theta_2)\, u(z)\, dz. \label{eqn:iae}\tag{1}\] Here, \(\phi_1\) is the encoding kernel, \(\phi_0\) a latent-space MLP, \(\phi_2\) the decoding kernel, and \(y \in \Omega_a\). The integrals are approximated as discrete sums, mapping variable-cardinality inputs to the fixed grid \(S_z\). This allows IAE-Net to natively handle different input resolutions. Blocks are connected with DenseNet-style skip connections.

The kernels \(\phi_1, \phi_2\) in equation 1 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 \(X_\text{bnd}\) while the decoder projects onto \(X_\text{dom}\), 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 \(X_\text{bnd} \in \mathbb{R}^{M \times d}\) be the boundary/sensor coordinates and \(X_\text{dom} \in \mathbb{R}^{N \times 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 \in \mathbb{R}^{K \times d_t}\); \(K\) is the total token count, \(d_t\) the token dimension.

Boundary feature tokenization. The boundary coordinates \(X_\text{bnd}\) are mapped to \(K_\text{bnd}\) tokens \(T_\text{bnd} \in \mathbb{R}^{K_\text{bnd} \times d_t}\) using a PointNet-style encoder [44], which applies shared MLPs followed by permutation-invariant global max-pooling, then projects and reshapes the resulting vector into \(K_\text{bnd}\) tokens: \[T_\text{bnd} = \text{Proj}\!\left(\max_{p \in X_\text{bnd}} \text{MLP}(p)\right).\]

State slice tokenization. Adapting the soft-clustering aggregation of Transolver [9] to inputs that include both spatial coordinates and physical values, we let \(X \in \mathbb{R}^{P \times d}\) be a point cloud (instantiating as \(X_\text{dom}\) for forward problems or \(X_\text{bnd}\) for inverse problems) and \(U\) the associated values. The concatenated state \([X, U]\) is mapped to \(K_\text{slice}\) soft-cluster tokens via: \[W = \text{Softmax}(\text{MLP}_\text{cluster}([X, U])), \quad H = \text{MLP}_\text{feat}([X, U]),\] \[T_\text{slice} = \text{diag}(W^\top \mathbf{1})^{-1} W^\top H,\] with mass-normalized aggregation; \(\mathbf{1} \in \mathbb{R}^P\) is a vector of ones. A linear projection brings \(T_\text{slice}\) to the standard token dimension \(d_t\). The initial token set \(T^{(0)} = [T_\text{bnd}, T_\text{slice}]\) 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. 1 ) so that the learnable kernel is conditioned on the geometry tokens \(T^{(l)}\). Spatial coordinates are first lifted via Fourier feature encodings \(\gamma(x) \in \mathbb{R}^{2 d_F + 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: \[C(x) = \text{Softmax}\!\left(\frac{(\gamma(x) W_Q)(T^{(l)} W_K)^\top}{\sqrt{d_k}}\right)(T^{(l)} W_V),\] where \(W_Q, W_K, W_V\) are learnable projections and \(d_k\) the scaling factor. The kernel is parameterized by an MLP taking the concatenated spatial features and geometric context: \[\label{eqn:cond} K(x, z; \theta) = \text{MLP}_\text{kernel}\!\left([\gamma(x), \gamma(z), U(x), C(x)]\right),\tag{2}\] where \(U(x)\) is the input feature at point \(x\). The integrals in (1 ) with this kernel are approximated as: \[I(z) = \frac{1}{|X|}\sum_{x_j \in X} K(x_j, z; \theta)\, U(x_j),\] followed by a learned residual connection and GELU non-linearity.

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 1 ) of IAE-Net, but with kernels conditioned on geometry tokens via cross-attention (Equation 2 ). 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(Nm)\) with \(m \ll 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 \in \mathbb{R}^{N \times c}\) is defined on \(X_\text{dom}\). The network applies \(L\) Adaptive IAE blocks with dense skip connections: the input \(H^{(l)}_\text{in}\) 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 \in \mathbb{R}^{M \times c_\text{obs}}\) at boundary locations \(X_\text{bnd}\), and the target lives on the dense interior mesh \(X_\text{dom}\). To bridge this gap, an Observation-to-Domain Decoder first transfers information from the boundary to a fixed set of \(K_\text{modes}\) learnable latent coordinates \(Z \in \mathbb{R}^{K_\text{modes} \times d}\), then projects from the latent set onto the interior mesh. Concretely, the boundary data is embedded as \(H_\text{bnd} = \text{MLP}_\text{in}(U) \in \mathbb{R}^{M \times w}\). A geometry-conditioned boundary integral lifts it to the latent set: \[H_\text{mode}(z_k) = \frac{1}{|X|} \sum_{x_j \in X_\text{bnd}} K_\text{enc}(x_j, z_k; T^{(0)})\, H_\text{bnd}(x_j)\, \qquad \forall z_k \in Z.\] A fixed-size MLP processes the latent representation, \(\tilde{H}_\text{mode} = \text{MLP}(H_\text{mode})\), after which a second geometry-conditioned integral projects onto \(X_\text{dom}\): \[H^{(0)}(y_i) = \frac{1}{|Z|}\sum_{z_k \in Z} K_\text{dec}(y_i, z_k; T^{(1)})\, \tilde{H}_\text{mode}(z_k)\, \Delta z_k, \qquad \forall y_i \in X_\text{dom}.\] The resulting \(H^{(0)} \in \mathbb{R}^{N \times 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 \(L^2\) error, \[\mathcal{L}(\theta) = \mathbb{E}_{f \sim \mu}\!\left[\frac{\|(\Psi_n(f; \theta) - \Psi(f)) \odot M\|_2}{\|\Psi(f) \odot M\|_2 + \epsilon}\right], \label{eq:masked-loss}\tag{3}\] where \(\odot\) 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 7.

4.1 Benchmarks↩︎

4.1.1 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 governed by \(\nabla \cdot (a({x}) \nabla u({x})) = 0\) on \(\Omega \subseteq \mathbb{R}^2\). 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 \(\Lambda_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 \(\sigma_s({x})\) of a tissue-like medium from boundary light measurements. The transport of photons through the medium \(\Omega \subset \mathbb{R}^2\) is governed by the stationary radiative transport equation (RTE): \[v \cdot \nabla_x u(x, v) + \sigma_t(x)\,u(x, v) = \sigma_s(x) \int_{S^{d-1}} \Phi(v \cdot v')\,u(x, v')\,dv', \quad x \in \Omega,\] where \(u(x,v)\) is the photon density at position \(x\) traveling in direction \(v\), \(\sigma_t = \sigma_s + \sigma_a\) is the total attenuation, and \(\Phi\) is the scattering phase function. The boundary observable is the Albedo operator \(\Lambda : L^1(\partial\Omega) \to L^1(\partial\Omega)\), which maps incoming illumination \(\phi(x)\) to the outgoing flux \(J_+(x) = \int_{v \cdot n_x > 0}(v \cdot n_x)\,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 \(\Lambda_0\) of a homogeneous baseline (\(\sigma_{s,0}{=}1\)) is used as a reference, and the inverse operator maps the difference \(\Delta\Lambda = \Lambda - \Lambda_0\) to the interior scattering perturbation \(\delta\sigma_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_\infty=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 \(-\Delta 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 \(\kappa( {x})\) to the pressure \(u( {x})\) satisfying \(-\nabla \cdot (\kappa \nabla 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 6.

4.2 Discussion↩︎

4.2.1 Inverse Problems↩︎

Table  1 reports median relative \(L^2\) 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.% for LNO), 21% on OT (1.41% vs.% for GINO), and 64% on the Airfoil task (0.58% vs.% 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.

4.2.2 Boundary Value Problems↩︎

Figure 4 visualizes representative reconstructions on some MCB shapes. On the 3D Poisson BVP (Table ¿tbl:tab:BVP?), GAIA outperforms all baselines on every shape category, with the largest gains on gears (3.44% vs.% for NGF) and screws & bolts (3.68% vs.% 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 2 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.% for Transolver). On 3D Darcy flow, Transolver achieves lower error than GAIA at the training resolution (0.73% vs.%). 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 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.

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 5 (a) 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 9.2, Figure 14).

4.4 Efficiency analysis↩︎

Table 3 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 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\times\) GAIA’s latency at \(40\times\) 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 \(\eta \in \{0\%, 1\%, 2\%, 5\%, 10\%\}\) using an additive Gaussian noise model scaled by the per-sample \(L^\infty\) norm. Figure 5 (b) shows graceful degradation across all three inverse benchmarks: at \(\eta = 10\%\), error grows by roughly \(13\times\) on Airfoil, \(4\times\) on EIT, and \(7\times\) 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 ¿tbl:tab:ablations? 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.% 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 (9-10) in Appendix 10.

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 discretization invariance. Thirdly, the geometry-conditioning cross-attention scales as \(\mathcal{O}(NK)\) in point-cloud size \(N\) and token count \(K\); for very large 3D point clouds (\(N > 10^6\)), 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.

6 Datasets and Benchmarks↩︎

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

6.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: \[\begin{align} &\nabla \cdot (a(x) \nabla u(x)) = 0, \text{ for } x \in \Omega \subseteq \mathbb{R}^d,\;d= 2,3, \label{eqn:laplace}. \end{align}\tag{4}\] Here, \(a(x)\) is the conductivity distribution of the medium. The goal is to recover this distribution using only the boundary data on \(\partial \Omega\) using the Dirichlet-to-Neumann map (DtN map, also known as the Calderon operator [46]) \(\Lambda_a\) given by: \[\Lambda_a: H^{-1/2}(\partial \Omega) \ni a \nabla u_{|\partial \Omega} \cdot n \mapsto u_{|\partial \Omega} \in H^{1/2}(\partial \Omega).\] where \(H^{-1/2}(\partial \Omega)\) is the space of bounded linear functionals on the Sobolev space \(H^{1/2}(\partial \Omega)\). In the case of full boundary measurement, it is a known result [46] that the inverse map \(\mathcal{F}^{-1}: \Lambda_a \mapsto 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 4 on strictly star-shaped domains \(\Omega\) defined by the polar radius \(r(\theta) = r_0[1 + c_1\cos(4\theta) + c_2\cos(8\theta)]\). 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: \[a(x, y) = \exp\left(\sum_{k=1}^{m} c_k \sin(k\pi x) \sin(k\pi y)\right),\] where \(m\) is drawn uniformly from \(\{1, \dots, 5\}\), and the coefficients \(\{c_k\}_{k=1}^m\) are sampled independently from a uniform distribution \(\mathcal{U}([-1, 1])\). For each sampled conductivity, the governing equation was solved numerically using quadratic Finite Element Method (FEM) under \(L=20\) distinct Dirichlet boundary conditions. The boundary data \(\{g_\ell\}_{\ell=1}^{L}\) corresponds to plane waves incident from equispaced angles \(\theta_\ell = \frac{2\pi (\ell-1)}{L}\), formulated as: \[u(x,y)\big|_{\partial \Omega} = g_\ell(x,y) = \cos\left(2\pi(x \cos \theta_\ell + y \sin \theta_\ell)\right).\] The corresponding Neumann measurement, \(\Psi_\ell = a \frac{\partial u}{\partial \nu}\), representing the current flux, was evaluated at \(M=272\) boundary sensors for each excitation \(\ell\). Figure 6 shows the Dirichlet and Neumann measurements on the boundary for a given star-shaped domain and its corresponding recovered internal conductivity distribution.

6.2 Inverse Problem: Optical Tomography↩︎

In optical tomography, the transport of photons through a scattering and absorbing medium occupying a bounded domain \(\Omega \subset \mathbb{R}^d\) is modeled by the stationary radiative transport equation (RTE). Let \(u(x, v)\) denote the particle density at spatial location \(x \in \Omega\) and velocity direction \(v \in S^{d-1}\). The medium is characterized by its scattering coefficient \(\sigma_s(x) \ge 0\) and absorption coefficient \(\sigma_a(x) \ge 0\). The RTE is given by: \[\begin{align} v \cdot \nabla_x u(x, v) + \sigma_t(x)u(x, v) &= \sigma_s(x) \int_{S^{d-1}} \Phi(v \cdot v')u(x, v')dv', \quad x \in \Omega, \\ u(x, v) &= \phi(x, v), \quad (x, v) \in \Gamma_-. \end{align}\] Here, \(\sigma_t(x) = \sigma_a(x) + \sigma_s(x)\) is the total attenuation coefficient, and \(\Phi(v \cdot v')\) is the scattering phase function, normalized such that \(\int_{S^{d-1}} \Phi(v \cdot v') dv' = 1\). The boundary phase space is partitioned into inflow (\(\Gamma_-\)) and outflow (\(\Gamma_+\)) boundaries, defined as \[\Gamma_\pm = \{(x, v) \in \partial \Omega \times S^{d-1} : \pm n_x \cdot v > 0\},\] where \(n_x\) is the unit outer normal vector at \(x \in \partial \Omega\). 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, \(\phi(x) \in L^1(\partial \Omega)\). Similarly, the measurements on the outflow boundary record the outgoing angular flux (or current), integrating the particle density over all outward-pointing velocities: \(J_+(x) = \int_{v \cdot n_x > 0} (v \cdot n_x) u(x, v) dv\). The associated boundary observation operator, known as the Albedo operator \(\Lambda\), thus maps the spatial input on the boundary to the spatial output: \[\Lambda : L^1(\partial \Omega) \to L^1(\partial \Omega), \quad \Lambda : \phi(x) \mapsto J_+(x).\] The continuous inverse problem seeks to recover the unknown medium properties, characterized by the scattering and absorption coefficients, from the knowledge of \(\Lambda\). Formally, this defines the continuous inverse map \(\mathcal{F}^{-1}\) from the space of bounded linear operators to the space of continuous medium parameters:\[\mathcal{F}^{-1} : \mathcal{L}\left(L^1(\partial \Omega), L^1(\partial \Omega)\right) \to C(\Omega) \times C(\Omega), \quad \mathcal{F}^{-1}(\Lambda) = (\sigma_s, \sigma_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 \(\Omega \subset \mathbb{R}^2\) to randomly scaled convex pentagons. We assume the absorption coefficient is a known, spatially invariant constant, \(\sigma_a(x) \equiv 0.01\), reducing the inverse problem to the recovery of the spatially varying scattering coefficient \(\sigma_s(x)\). The highly forward-peaked anisotropic scattering in biological tissues is modeled using the 2D Henyey-Greenstein (HG) phase function,\[\Phi(v \cdot v') = \frac{1 - g^2}{2\pi(1 + g^2 - 2g(v \cdot 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: \[\sigma_s(x) = 1.0 + \left({1 + \exp\left(-3 \sum_{k=1}^m c_k \sin(\omega_{x,k} \pi x'_1) \sin(\omega_{y,k} \pi x'_2)\right)}\right)^{-1},\] where the number of modes \(m\) is drawn uniformly from \(\{2, \dots, 5\}\), the expansion coefficients \(c_k\) are sampled independently from a uniform distribution \(\mathcal{U}([-1, 1])\), and the spatial frequencies \(\omega_{x,k}\) and \(\omega_{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 \(\sigma_s(x) \in (1.0, 2.0)\).

The continuous boundary illumination is discretized using \(N_{src} = 32\) localized isotropic point sources distributed uniformly along \(\partial \Omega\). The forward transport problem is solved using the Discrete Ordinates (\(S_{32}\)) 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 \(N_{rec} = 32\) interleaved receiver locations. This discrete setup approximates the Albedo operator as a dense matrix \(\Lambda \in \mathbb{R}^{32 \times 32}\). To isolate the scattering anomalies and suppress structural geometric artifacts, we employ a difference imaging approach. A baseline Albedo matrix \(\Lambda_0\) is simulated using a known, constant background scattering coefficient \(\sigma_{s,0} = 1.0\). The objective of the discrete inverse problem is then to infer the scattering perturbation \(\delta \sigma_s(x) = \sigma_s(x) - \sigma_{s,0}\) strictly from the measured difference data \(\Delta \Lambda = \Lambda - \Lambda_0\).

6.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, \[\frac{\partial \rho_f}{\partial t} + \nabla \cdot (\rho_f {v}) = 0, \qquad \frac{\partial \rho_f {v}}{\partial t} + \nabla \cdot (\rho_f {v} \otimes {v} + p {I}) = 0, \qquad \frac{\partial E}{\partial t} + \nabla \cdot \bigl((E + p) {v}\bigr) = 0,\] where \(\rho_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 \(\rho_\infty = 1\), \(p_\infty = 1.0\), \(M_\infty = 0.8\), \({AoA} = 0^\circ\), with a no-penetration condition imposed at the airfoil surface. Here \(M_\infty\) 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 \(N_x \times N_y = 221 \times 51 = 11{,}271\) quadrilateral elements, refined near the airfoil surface. Airfoil geometries are parameterized via a design element approach [53] in which the baseline NACA-0012 profile is deformed by displacing the control nodes of an enclosing cubic design element in the vertical direction, with displacements drawn from \(d \sim \mathcal{U}[-0.05, 0.05]\). The dataset is partitioned into \(N_{\mathrm{tr}} = 1{,}000\) training, \(N_{\mathrm{val}} = 100\) validation, and \(N_{\mathrm{te}} = 100\) test configurations.

We formulate it as a sparse-to-full field reconstruction problem. The quantity of interest is the steady-state Mach number field \(u : \Omega \to \mathbb{R}\), discretized over all \(11{,}271\) mesh nodes. At inference time, only a sparse subset of observations \(\{u( {x}_j)\}_{j \in \mathcal{O}}\), with \(\mathcal{O} \subset \{1, \ldots, N_x N_y\}\), is available as input, and these observations are further corrupted by zero-mean Gaussian noise scaled to the per-sample peak Mach number: \[\tilde{u}_i( {x}_j) \;=\; u_i( {x}_j) \;+\; \eta_{i,j}, \qquad \eta_{i,j} \sim \mathcal{N}\!\left(0,\, \sigma^2 \|u_i\|_\infty^2\right), \quad j \in \mathcal{O}, \label{eq:obs-noise}\tag{5}\] with \(\sigma = 0.01\), corresponding to \(1\%\) of the per-sample \(L^\infty\) 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.

6.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 \(f_{\text{analytical}}(x,y,z)\). The base function is constructed using trigonometric variations along all three spatial axes:\[f_{\text{analytical}} = T_1(x,y) + T_2(x,y) + T_3(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:\[T_1(x,y) = -\left((A\pi)^2 + (AC\pi)^2\right) \sin(A\pi x) \cos(AC\pi y)\]\[T_2(x,y) = (A\pi)^2 \cos(A\pi x)\left(1 - \sin(AB\pi y)\right) + (AB\pi)^2 \sin(AB\pi y)\left(1 - \cos(A\pi x)\right)\]\[T_3(z) = 2(AD\pi)^2 \cos(2AD\pi 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(x^3 - 3xy^2) + F(y^3 - 3x^2y) + (x^2 - z^2)\]For each generated sample, the coefficients \(E\) and \(F\) are sampled from uniform distributions such that \(E \sim \mathcal{U}(-1.0, 1.0)\) and \(F \sim \mathcal{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.

6.5 Forward Problem: 3D Darcy Flow↩︎

We consider the stationary Darcy flow through a heterogeneous porous medium. Let the spatial domain \(\Omega \subset \mathbb{R}^3\) represent three-dimensional domain. The fluid pressure \(u( {x})\) is modeled by the elliptic partial differential equation: \[-\nabla \cdot (\kappa( {x}) \nabla u( {x})) = 1, \quad {x} \in \Omega,\] subject to homogeneous Dirichlet boundary conditions, \(u|_{\partial \Omega} = 0\). Here, \(\kappa( {x})\) represents the strictly positive permeability field of the medium. The operator learning objective is to approximate the continuous solution operator \(\mathcal{G}_{d}\) that maps the permeability field, defined on a variable geometries, to the resulting scalar pressure field. Formally, for a given domain \(\Omega\), the mapping is defined as: \[\mathcal{G}_{df} : L^\infty(\Omega) \to H_0^1(\Omega), \quad \mathcal{G}_{df}(\kappa) = u.\]

To learn this operator, we construct a dataset by solving the Darcy flow equation on \(\Omega \subset \mathbb{R}^3\). The boundary \(\partial \Omega\) of each sample is generated as a "star shape" defined by a random spherical harmonic expansion up to degree \(l_{max}=4\): \[r(\theta, \phi) = \sum_{l=0}^{4} \sum_{m=-l}^{l} c_{lm} Y_l^m(\theta, \phi), \quad c_{lm} \sim \mathcal{N}(0, \sigma_l).\] The resulting geometries are discretized into watertight unstructured tetrahedral meshes. The permeability field \(\kappa( {x})\) is computed as a continuous log-normal random field, \(\kappa( {x}) = \exp(\gamma \cdot \mathcal{S}( {x}) / \sigma_{\mathcal{S}})\), where \(\mathcal{S}( {x})\) is a superposition of \(M=30\) random Fourier modes: \[\mathcal{S}( {x}) = \sum_{j=1}^{M} \alpha_j \cos( {k}_j \cdot {x} + \phi_j).\] Here, the wave vectors \({k}_j \sim \mathcal{N}( {0}, \sigma_k^2 {I})\) control the spatial correlation length, phases \(\phi_j \sim \mathcal{U}[0, 2\pi]\) ensure translation invariance, and amplitudes \(\alpha_j \sim \mathcal{N}(0, 1)\) provide variance. The field is normalized by its sample standard deviation \(\sigma_{\mathcal{S}}\) and scaled by \(\gamma=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 (\(P_1\)) basis functions for both the solution space and the coefficient field. Following the FEM solution, we extract the training data by randomly sampling \(N_p=7000\) points within the domain volume for each simulation.

6.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: \[\rho \partial_{tt} {u} - \nabla \cdot \boldsymbol{\sigma} = {0},\] where \(\rho\) represents the mass density, \({u}\) the displacement field, and \({\sigma}\) 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 \(\sigma\), which is available at \(972\) unstructured coordinates. The operator learning task thus seeks to approximate the continuous solution operator \(\mathcal{G}_{el}\) 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) = \text{dist}(x, \Gamma_{\text{inner}})\). The formal operator mapping is thus defined as:\[\mathcal{G}_{el} : C(\Omega) \ni d \to \boldsymbol{\sigma} \in L^2(\Omega; \mathbb{R}^{2 \times 2}).\]

Figure 11 (a) 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.

6.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 \(\Omega = (0,1)^2\). The governing PDE is given by: \[\Delta u(x) = f(x),\;x \in (0,1)^2 ,\] subject to homogeneous Dirichlet boundary conditions, \(u|_{\partial \Omega} = 0\). The objective of the neural operator is to learn the inverse Laplacian operator \(\mathcal{G}_{pg} = \Delta^{-1}\), which maps the right-hand side source function to the scalar solution field. Formally, this defines the continuous solution operator: \[\mathcal{G}_{pg} : L^2(\Omega) \to H_0^1(\Omega) \cap H^2(\Omega), \quad \mathcal{G}_{pg}(f) = u.\] The source term \(f\) is procedurally generated as a superposition of \(M\) random Gaussian pulses: \[f(x, y) = \sum_{i=1}^{M} \exp\left( -\frac{(x - \mu_{x,i})^2 + (y - \mu_{y,i})^2}{2\sigma_i^2} \right),\] where \(M\) is an integer drawn from a geometric distribution. The mean coordinates are sampled uniformly such that \(\mu_{x,i}, \mu_{y,i} \sim \mathcal{U}(0,1)\), and the standard deviations are sampled as \(\sigma_i \sim \mathcal{U}(0.025, 0.1)\). The ground truth data pairs \((f, u)\) are evaluated on a high-resolution \(128 \times 128\) uniform grid. We refer to [36] for further dataset details. Figure 11 (b) visualizes a representative input source function and its corresponding computed solution.

7 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 \(\omega_0 = 30\). All hyperparameters for experiments are given in Table 5.

3pt

8 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.

8.1 CORAL↩︎

CORAL [11] casts operator learning as a regression in the latent spaces of two implicit neural representations (INRs). Each input sample \(a_i\) and output sample \(u_i\) is encoded by fitting a shift-modulated SIREN via auto-decoding: a shared base network is trained jointly with per-sample latent codes \(z_{a_i}, z_{u_i}\), 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 \(z_{a_i} \mapsto z_{u_i}\) in latent space. Because the INRs are queried point-wise, CORAL natively supports varying input and output meshes across samples.

We use the official implementation³. 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 \(z_{a} \mapsto z_{u}\) is a 3-block residual MLP with Swish activations, trained with Adam (no weight decay) and a step decay of \(\gamma=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.

8.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 implementation⁴. 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 6.

8.3 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 \(\mathcal{O}(N^2)\) to \(\mathcal{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 implementation⁵. 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.

8.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 ⁶. For 2D problems, the latent token 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.

8.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 repository⁷ 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.

8.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 repository⁸ without modification.

9 Additional Results↩︎

9.1 Noise Robustness↩︎

We adopt a relative-Gaussian noise model scaled by the per-sample \(L^\infty\) norm of the measurement field. Let \(x\) denote a clean measurement vector; the corrupted measurement is: \[\label{eq:noise95model} \tilde{x} = x + \varepsilon, \qquad \varepsilon \sim \mathcal{N}\!\left(0,\; \bigl(\eta\,\|x\|_\infty\bigr)^2 I\right).\tag{6}\] for noise levels \(\eta \in \{0\%, 1\%, 2\%, 5\%, 10\%\}\)

The three benchmarks differ in what the measurement vector \(x\) is and at what granularity the \(L^\infty\) norm is computed.

• Airfoil. The measurement is the partial Mach field observed at a fraction of mesh points. A single \(\|x\|_\infty\) 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 \(\Delta\Lambda \in \mathbb{R}^{N_{rec} \times N_{src}}\). Because receivers have very different outgoing flux magnitudes, we use a per-receiver-row \(L^\infty\) scale: \(\varepsilon_{r,s} \sim \mathcal{N}(0, (\eta\,\|\Delta\Lambda_{r,:}\|_\infty)^2)\).

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 \(L^2\) error across these seeds; results are summarized in Table 7. 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.

9.2 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.

9.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 \times 128\) grid and subsampled to \(64 \times 64\) for standard training. For augmented training, the \(64 \times 64\) data is upsampled back to the canonical \(128 \times 128\) resolution using bicubic interpolation. Multi-resolution views are then generated by striding this grid. Each training iteration simultaneously optimizes over multiple resolutions \(\mathcal{R} = \{32 \times 32, 32 \times 64, 64 \times 32, 64 \times 64, 64 \times 128, 128 \times 64, 128 \times 128\}\), with the total loss: \[\mathcal{L}_{\text{total}} = \sum_{T \in \mathcal{R}} \lambda_T \cdot \mathcal{L}(\hat{u}_T, u_T),\] where \(\lambda_T\) are per-resolution weights (set to 1.0 in this experiment). As shown in Figure 15, multi-resolution 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.

9.4 Efficiency analysis↩︎

Table 3 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.

9.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 8 reports the resulting statistics. The small standard deviation indicates stable performance across runs. For airfoil, we remove the random noise to isolate model randomness.

10 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.

10.1 Token Dimension↩︎

Table 9 reports the median relative \(L^2\) 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.

10.2 Number of Geometric Tokens↩︎

Table 10 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.

10.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 11 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.

10.4 Number of Blocks↩︎

Table 12 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.

10.5 Geometry Conditioning Ablation Details↩︎

Table ¿tbl:tab:ablations? 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.

10.6 Tokenizer Ablation Details↩︎

The IAE refinement stage stacks several adaptive blocks, each applying integral transforms followed by token-updated cross-attention. Table ¿tbl:tab:ablations? evaluates the contribution of each token pathway by removing either the boundary feature tokens (produced by the PointNet-style encoder from \(X_{\text{bnd}}\)) 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.

References↩︎

---

1. Significant technical contribution.↩︎

2. 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.
   Correspondence: Meenakshi Krishnan (mkrishn9@umd.edu).↩︎

3. https://github.com/LouisSerrano/coral↩︎

4. https://github.com/neuraloperator/neuraloperator↩︎

5. https://github.com/thuml/Transolver↩︎

6. https://github.com/camlab-ethz/GAOT↩︎

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

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