RIGNO: A Graph-based Framework For Robust And Accurate Operator Learning For PDEs On Arbitrary Domains
Authors: Sepehr Mousavi
Shizheng Wen
Levi Lingsch
Maximilian Herde
Bogdan Raonic
Siddhartha Mishra
Venue: NeurIPS 2025 poster ()
Keywords: scientific machine learning
operator learning
partial differential equation
neural operator
graph neural network
AI for science
computational science
Primary area: machine_learning_for_sciences
TL;DR: We propose a very performant graph-based neural operator architecture for learning the solution operator of PDEs from data on arbitrary domain discretizations, which has been tested on a challenging suite of benchmark datasets.

Abstract: Learning the solution operators of PDEs on arbitrary domains is challenging due to the diversity of possible domain shapes, in addition to the often intricate underlying physics. We propose an end-to-end graph neural network (GNN) based neural operator to learn PDE solution operators from data on point clouds in arbitrary domains. Our multi-scale model maps data between input/output point clouds by passing it through a downsampled regional mesh. The approach includes novel elements aimed at ensuring spatio-temporal resolution invariance. Our model, termed RIGNO, is tested on a challenging suite of benchmarks composed of various time-dependent and steady PDEs defined on a diverse set of domains. We demonstrate that RIGNO is significantly more accurate than neural operator baselines and robustly generalizes to unseen resolutions both in space and in time. Our code is publicly available at github.com/camlab-ethz/rigno.

RIGNO: A Graph-based Framework For Robust And Accurate Operator Learning For PDEs On Arbitrary Domains
Sepehr Mousavi
Shizheng Wen
Levi Lingsch2.
Maximilian Herde
Bogdan Raonic
Siddhartha Mishra
Computational Mechanics Group, ETH Zurich, Switzerland
2 Seminar for Applied Mathematics, ETH Zurich, Switzerland ETH Al Center, Switzerland
Abstract
Learning the solution operators of PDEs on arbitrary domains is challenging due to the diversity of possible domain shapes, in addition the often intricate underlying physics. We propose an end-to-end graph neural network (GNN) based neural operator to learn PDE solution operators from data point clouds in arbitrary domains. Our multi-scale model maps data between input/output point clouds by passing it through downsampled regional mesh. The approach includes novel elements aimed at ensuring spatio-temporal resolution invariance. Our model, termed RIGNO, is tested various time-dependent and steady PDEs defined on diverse set of domains. We demonstrate that RIGNO is significantly more accurate than neural operator baselines and robustly generalizes to unseen resolutions both in space and in time. Our code is publicly available atgithub..o//mmmlldeezziiggg
challenging suite of benchmarks composed of
Introduction
Partial Differential Equations (PDEs) mathematically model wide variety of interesting phenomena physics and engineering, making their study a matter of great scientific importance (9]. Given the paucity of explicit solution formulas, numerical methods [36), such finite difference, finite element and spectral methods constitute the main tools in the quantitative study of PDEs. The computational cost of existing numerical methods or solvers for PDEs can be prohibitively expensive, particularly in the context of problems such as uncertainty quantification, inverse problems, control, and design, where PDE solvers have to be called multiple times 35 This necessitates the design of fast and accurate surrogates for PDE solvers.
In this context, machine learning (ML) algorithms are increasingly being used surrogates for PDE solvers, see [31 and the references therein. In particular, operator learning, i.e., where the task at hand is to learn the underlying solution operator of a PDE from data, has emerged as the dominant paradigm in the application ML techniques for PDEs. A large variety of operator learning algorithms have recently been proposed, with examples being DeepONets 7 28 Fourier Neural Operator (FNO) [20]. Low-rank Neural Operators [22]. Convolutional Neural Operators (CNO) (8] and transformer-basee
operator learning algorithms [6] 34 12 13 D. among many others. The
success of these approaches may largely be contributed to their multi-scale learning. For example, the FNO learns frequencies of different scales within the data, CNO learns by down- and up-sampling
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

via convolutions, and transformer-bee based algorithms often use multi scale windowing. As a result, these algorithms learn both short- and long-range physical dependencies of the solutions of PDEs.
An overwhelming majority of these operator learning algorithms or neural operators [66 assume that the underlying domain is Cartesian and is discretized on a regular grid; however, in practice, and particularly in engineering, PDEs are posed on domains with very complex and highly non- Cartesian geometries. Thus, the direct application of many of the aforementioned operator learning algorithms is difficult or impractical. Therefore, we require algorithms that can learn the solution operators of PDEs on arbitrary domains. Various approaches have been proposed to extend neural operators to arbitrary domains. To start, methods such as the FNO and the CNO, originally designed for
Cartesian domains, can be extended to non Can-Cartssian domains via domain masking [38 [3]. Another approach for extending FNO to arbitrary domains is to directly evaluate the discrete Fourier transform rather than FFT [44 Moreover, one can transform the underlying problem to a Cartesian domain
either through extensions [1], learned diffeomorphisms as in Geo-FNO [9], or through learned interpolations between non- non-Cartesian physical domain and a Cartesian latent domain via graph neural operators [21 as in GINO [3 and references therein. Applying attention mechanisms to downsampled representation, as demonstrated in UPT I and Transolver [45], offers another viable approach. A conceptually distinct line of work leverages implicit neural representations 4ll 42). which accommodate arbitrary geometries by encoding input and output functions into a shared latent vector that can be continuously queried at any spatial coordinate. The mapping between input and output functions is then learned in this latent space through a separate training stage.
An alternative approach to learn operators for PDEs on arbitrary domains is provided by graph neural networks (GNNs) [5] This is an ideal MI. framework for handling relational data such as molecules, networks, or unstructured computational grids. The use of graph neural networks for simulating PDEs has received widespread attention 44
8] to name a few. Particularly, approaches such as continue to demonstrate success in learning physical simulations. Graphs are able to implicitly capture the challenging collocation point distributions of unstructured meshes and arbitrary point clouds, making methods such as MeshGraphNet ideal for the non- Cartesian domain. However, genuine learning of operators, rather than a specific finite- dimensional discretization of them, requires the ability of the underlying algorithm to display some form of space-time resolution invariance [66] 2). i.c., the operator learning algorithm should be able to process inputs and outputs on different spatial and temporal resolutions, both during training and testing, without any decrease in its overall accuracy. Existing graph-basee approaches such as MeshGraphNet cannot satisfy this key property, and even many neural operators on Cartesian grids do not inherently capture temporal resolution invariance. Moreover, due to the local nature of their graph construction, MeshGraphNe struggle to capture long-range dependencies within the data. Therefore, given the success of multi-scale neural operators such as FNO and CNO, the question naturally arises: does there exist a complementary, multi-scale, graph-based based approach which is capable of learning space-time resolution invariant operators on complicated distributions of input points?
39
33
MeshGraphNet [33] and its derivatives 110 77
In this work, we answer this question affirmatively by proposing RIGNO, Region Interacting Graph Neural Operator. RIGNO is an end-to-end graph neural network based on the encode-process- decode paradigm 44 39. relying an algorithm to construct multi-scale regional mesh. We also incorporate several novel elements which, in combination with the algorithmic developments, push this algorithm achieve state-of-the-art performance. The contributions of this work are as follows:
The introduction of flexible regional mesh which simultancously incorporates resolution invariance and multi-scale feature learning.
An edge masking technique for GNNs in the context of scientific machine learning, which improves resolution invariance, facilitates optimization during training, and can provide model uncertainty estimates at inference.
The introduction of temporal fractional pairing, a fine-tuning strategy for time-continuous neural operators that unlocks generalization to higher time resolutions.
Extensive experimental results that demonstrate significant improvements to the state-of- the-art across variety of benchmarks with domain geometries including structured and unstructured meshes, random point clouds, and Cartesian grids.

---

2 Methods
Problem formulation. We consider the following generic time-dependent PDE with S2 = ()) the time domain, 22 C E the d-dimensional spatial domain, u c(I;; X) the solution of the PDE with X' C (Q2 R²) for some 1 S p < oc, c E Q spatially varying coefficient, F differential operator, B boundary operator, and E X' the initial condition.
auu F(c,t, u, Vau, V.n,...) Bu(t, x) = 0, u(0, x)) = a(x):
v(t, E Q x J2z, v(t,z) E lQ x 002, Vr E fΩg.
(I)
The solutions of the PDE 1 are given in terms of the solution operator S' X XXX + X' which maps the initial conditions and coefficients to the solution u(t) at any time. With the semi-group property for solutions of time-dependent PDEs [9], we can define an extended solution operator
G : X Q x S2 x R + X. which maps the solution u(t) at any time and the coefficient c the solution utt + T) at later time & + T through the action,
('((()),,447)) = S (u(), u(t + +T).
(2)
Our goal here is to learn this extended solution operator G from data. Note that plugging t = 0, T 1, u()) = a yields St, recovering the solution operator of the PDE T We would also like to point that the (extended) solution operator depends implicitly on the underlying spatial domain Ωz which in turn, can change in some operator learning tasks.
Model architecture. As shown in Figure | our model is based on the well-known encode-procee process- decode paradigm [40 39). All the three components, i.e., encoder, processor, and decoder are graph neural networks (GNNs). The methods presented in this section are implemented using JAX (3) and learn the (extended) solution operator G (2). the inputs to our model are, in general, i) a set of coordinates in the domain == {x, E }i where Gg is the spatial domain, ii) u == {u, E IR",, }., discretized values of a function at known coordinates x, ili ) discretized values of the spatial coefficient at x coordinates, iv) the time & and v) the lead time т. The model's 's output r == {r; E IR') } represent discretized values of a parameterized map at the coordinates and is obtained by
are available on Gut As the aim
e := {c:
c,,,
Underlying graphs. The map To is a GNN and will process data that are defined on graphs. As shown in Figure (for an idcalized wwo-dimensional domain) and following ideas in [7]. we consider two underlying sets of nodes. The first set of nodes corresponds to the unstructured point cloud on which the input (and output) data is defined. We call this point cloud the physical mesh and the nodes that constitute as physical nodes. We do not perform any explicit information processing on these physical nodes. Most of our processing will take place in another bespoken set of nodes that we call as regional nodes and the point cloud constituted by them as the regional mesh. Here, the regional nodes are constructed by randomly sub-sampling the set of physical nodes with a sampling ratio which is a hyperparameter of our model. A graph is constructed to connect the regional nodes (see Figure T via bidirectional regional mesh edges. Moreover, we need to pass information between the physical nodes and the regional nodes, necessitating the construction of graphs between them. To this end, we connect physical nodes regional nodes (Figure T by unidirectional edges. These edges are constructed such that each regional node constitutes circular nearby physical nodes and serves to accumulate information from them. Based on the local density of regional nodes, minimum support radius is calculated for each of them such that their union covers the whole domain. Information from the regional mesh can be transmitted to the physical mesh in a similar way by constructing a graph with unidirectional edges. As shown in Figure [ these edges are constructed to accumulate information from multiple regional nodes into one physical node. Thus, construction of these unidirectional edges provides a very flexible upscaling and downscaling strategy for our model.
containing subset of
Needless to say. when possible, the graph construction must be adapted to the boundary conditions of the PDE I Once the underlying graphs have been constructed, we must also provide geometric information to the GNNs. We normalize the given coordinates and prepare them for all nodes and all
3

---

Edge masking. Considering the graph structure the encoder as an example, each regional node aggregates information from multiple physical nodes. For high-resolution inputs, many of these physical nodes can be redundant, as the underlying function can often be accurately reconstructed (c.g., through interpolation) using a smaller subset of nodes. During training, the model may implicitly disregard some of these redundant nodes and instead focus on a subset of dominant or central connections. However, at test time, when the input resolution differs, these dominant nodes may no longer exist, resulting in degraded performance. We identify this phenomenon as primary cause of the resolution sensitivity commonly observed in GNN-based operator learning models. To address this, we use random masking of edges in all message-passinng blocks of RIGNO. Before each message passing step (e.g., step p of the processor), each edge is masked with given probability block. In the next message passing step Feed forward
Poscce exc
MLP
Bradarelg c
Tayerhorn
Deaddrra elgi
FAC
Mewage pasiing
Message pooding
MS
Figure 1: General schematic of the RIGNO architecture for idcalized two-dimensional domain. The inputs are first independently projected to a latent space by feed-forward blocks. The information on the original discretization (physical nodes) is then locally aggregated to coarser discretization (regional nodes). Regional nodes are connected to each other by edges with multiple length scales. Several prsssage-pssiing steps are then applied on the regional nodes which constitute the processor. The processed features are then transmitted back to the original discretization by using similar edges in the encoder, before being independently projected back the desired output dimension via feed-forward block without normalization layers.
edges accordingly. We postpone the rather technical description of graph construction, geometric features, and incorporation of boundary conditions, particularly the novel additions to the model design to enforce periodic boundary conditions, to SM Section A Once the inputs, graphs and geometric information are available, the first step is to lift the initial features of all nodes and all edges to suitable high-dimensional latent spaces. To this end, feed-forward blocks are applied independently on the initial features each node and edge. Each feed-forward block consists of a shallow multi-layer perceptron (MLP) followed by normalization layers (see Figure [ Next, model processes this high-dimensional latent information by employing the following constituent GNNs.
Main components. The encoder is designed to transmit information from the underlying physical mesh (representing discretization of the domain Ω2)) to the coarser regional mesh. This is done by performing single message passing step with residual updates. The processor then operates on the latent information defined over the regional mesh by performing P sequential message passing steps. We incorporate residual updates and skip connections to maintain gradient flow and promote stable training for large P. The decoder serves to transfer processed information from the regional mesh back to the physical mesh. The first step in the decoder a single message-passing block. In the second step of the decoder, the physical node feature vectors independently scaled down the desired dimension to form the output r. Summarizing, as our model involves the learned interaction between nodes that each represent aggregated information from an underlying sub-region, call as Region Interaction Graph Neural Operator or RIGNO. More details about the GNN components are provided SMB
and does not participate in that specific

---

(e.g., step p + 1 of the processor), a different set of edges is randomly masked. Through random masking of the edges during training, the model cannot consistently rely on the same high-correlation connections. In the absence of these edges, the model is forced to use the weaker correlations. Additionally, learning from weaker or secondary correlations the processor enables the model to better integrate information across multiple spatial scales, rather than overfitting to dominant patterns in the training resolution. Although loosely inspired by dropout [43] edge masking is conceptually very different from it. While dropout masks trainable parameters of MLPs, edge masking disables pairwise interaction of spatial entities (nodes). This allows the model to discover weaker spatial correlations in the data. The set of the masked edges randomly changes in every epoch, which allows a single message- passgenga block to learn from all the edges during training. Since the messages computed from the neighbors are averaged, there is need for any further configuration training or inference. It is also possible to use a different masking probability or to completely disable it at any point during training or inference. As shown in the next section, edge masking enables us i) make RIGNO (approximately) resolution invariant, ii) obtain significantly lower errors, iii) reduce computational complexity by sparsely sampling edges during each step, and iv) provide a method for quantifying uncertainty in the predictions of RIGNO.
Temporal processing. We recall that our goal is to approximate the extended solution operator (u(()),,t,, (2), which is function of the input time and the lead time Both of these scalars are provided to the encoder of RIGNO as inputs and embedded into suitable high-dimensional latent spaces. However, merely providing lead time as an input may not be sufficient to learn the solution operator in deep models. Following [13. we introduce lead time as an additional conditioning in each hidden layer of RIGNO (see Figure D. We do this by adapting the normalization layers employed within our model. Let be an output vector for any layer and normalize it by an appropriate normalization layer (e.g.. BatchNorm, LayerNorm) to obtain the normalized output y. We further condition y with the lead time r by setting y = (1 (T(()) O y + TX(T), with À being small learnalle MLPs.
Finally, we need to provide a time marching strategy to update our model in time. In this context, RIGNO can be employed in very general time marching mode by setting,
Go dt,,, T) = auu + (To .'',t, ,t,,,),
where a, S are scalar parameters. To is the output of RIGNO, and Go 2 G is the final approximation to the extended solution operator (2) of the PDE I Setting a = 0 and 3 1 directly recovers an approximation to the output of the extended solution operator. In this case, RIGNO can be interpreted to directly provide the solution of the PDE ① at time + T. given the solution at time t, coefficient e and the lead time T as inputs. Similarly, setting a = 3 = 1 yields the residual of the solution at the later time, given the solution at current time. On the other hand, setting a = 1 and / the model output To approximating the time derivative of the solution operator as 3 is simply an approximation of the Taylor expansion in time of the solution operator, truncated at leading order. A clear advantage of this time-derivative update is that it automatically enforces consistency for r 0. We will explore these choices for time updates in our numerical experiments to determine which ones are the most suitable.
leads to
Loss functions and training strategy. For simplicity of exposition, here, we fix final time T and consider a fixed discretization Q2 = U,= [I, 1, (n). We assume that training samples are provided as trajectories (wm (t,))) for 1 < m < M, with each thm being the solution to the PDE ! with initial condition am and coefficient Cm-- Given these training trajectories and the fact that our model can be evaluated at any time, there are multiple ways in which the RIGNO model, denoted here by Go can be trained. We follow the recently proposed all2all training strategy (13) as training outut- pairs scale quadratically with respect to the number of time steps for time-dependent problems. Hence, the goal of training is to find parameters 0 E θ such that the following loss function can be be minimized,
J = MN
Σ
C. (umm(tn) af (f.)), Tff (fr.) Go (Mon (ft)) Cmn be, to fx))
(4)
with N = (N + 1)(N 2) and L being a measure of the mismatch between the ground truth and the model prediction, for instance either the mean absolute or the mean squared error.
5

---

A problem with this approach is that the smallest prediction lead time Tiin in the training set is given by the minimum of f,, f,-1 for all n. Hence, at inference, predicting for lead times T Tiin could lead to distribution shifts as smaller lead times have not been seen during training. Yet, fundamental requirement for a time-dependent neural operator is to predict arbitrarily small times correctly in order to ensure time confinuity. We tackle this problem by proposing to further fine-tune a RIGNO model, trained with the loss function 4 by using a modified loss function where ih (tn) in 4 is obtained by starting from an intermediate input function: Go (Tmmmtt.) Cm. t,, in - t.). Thus, our aim is to use solutions at intermediate time steps, which are not available in the training set, to make predictions at later time steps, available in the training data. However, we must find an approximation (t.) Rs Um (t.) ) as um ) is not available during training. Our solution is to use the pretrained RIGNO itself by setting ū(t.) Go (u(tæ), C, tk. t. tx). This self-consistent approach is reasonable as the modified loss function is only used during fine-tuning. Recycling model outputs as inputs has been employed in previous implementations for the fixed-step setting [4 20] to improve rollout stability. Here, we extend this concept to the time-continuous setting, aiming to generalize to arbitrarily small time steps. As shown in SM Figure C.3 the flow of the gradient can be stopped at the intermediate time i, for training efficiency. The utility of this fractional pairing strategy will be tested subsequently.
Given that the output time t and lead time r are both inputs to the trained RIGNO, a variety of approaches can be used for temporal processing inference. As in [33, to obtain the solution ut)
at time t, from an initial datum (0) = a, we can either directly apply RIGNO with lead time T = or we can employ different autoregressive rollout strategies as described in the SM Section C.2
3 Results
In this section, we demonstrate the performance of RIGNO on a challenging suite of numerical experiments by comparing it with well-established baselines. We also highlight the factors that underpin the observed performance of RIGNO.
RIGNO accurately learns the solution operators of PDEs on arbitrary domains. To test the performance of RIGNO for PDEs on arbitrary domains, we choose two types of datasets. First, we consider 4 datasets on point clouds obtained from structured and unstructured meshes, namely Heat-L-Sines (heat equation on L-shaped domain), Wave-C-Sinel (wave equation, propagation of waves in homogeneous medium of circular shape), AF (compressible transonic flow past airfoils taken from 19), and Elasticity (hyper-elastic deformations of varying domain shapes, also taken from [19). Next, we test how well RIGNO learns operators defined on random point clouds. To set this up, we choose 9 datasets from the PDEGym database [I3, which deal with four different PDEs; the incompressible Navier-Stokes equations, the compressible Euler equations, the Allen-Cahn equation, the wave equation on a spatially varying medium, and the Poisson equation. The datasets all originally defined uniform Cartesian grids. To obtain random point cloud, randomly subsample 9216 points from the available coordinates on a 128x 128 uniform grid, see SM Figure D.I Our dataset collection spans 8 PDEs and comprises 3 time-independent and 10 time-dependent problems. See SM Section D for more details about the datasets and their abbreviations used here.
In T we present the test errors on these 13 unstructured datasets 32 with two sizes of RIGNO, namely RIGNO-12 12 and RIGNO-18, with 12 and 18 layers in their processor modules, respectively. To contextualize the results, we consider five baselines, namely MeshGraphNet 33 Geo-FNO [19 FNO DSE 224, GINO (23), and UPT 33 All of these baselines are capable of handling inputs and outputs on arbitrary point clouds, yet use different routes to achieve this. More details about the basclines, including the selection of hyperparamcters, are outlined in SM Section G We observe from Table [ that RIGNO-18 outperforms all the basclines in all 13 datasets. Likewise, RIGNO-12 is the second best- persoperfor model for 10 of the 13 datasets. For Wave-C-Sines, Poisson-Gauss, and Elasticity, other architectures claim second best; yet, RIGNO-12 remains comparable in performance to RIGNO-18 18. Moreover, in some cases, the RIGNO variants are almost an order of magnitude more accurate than the baselines, showing excellent performance on these benchmarks.
Often, models are designed to have flexible inputs/outputs by potentially reducing their performance on Cartesian domains. To further illustrate the expressive power of RIGNO, we consider 9 datasets from the PDEGym database on their original Cartesian grids. this case, we can compare RIGNO with state-of-te the neural operators such as FNO, CNO, and ScOT, which are tailored for Cartesian

---

Table I: Benchmarks datasets with unstructured and Cartesian-grid space discretizations. Lowest and second lowest errors are highlighted. The errors for time-depenn siee-dependent problems are the lowest among three autoregressive strategies (see SM Section C.2 All training datasets have 1024 trajecto- ries/samples except for the AF dataset (2048). MGN corresponds to the MeshGraphNNr baseline.
Dataset
Median relative L' error [%]
(point cloud)
RIGNO-18
RIGNO-12
MGN
Geo-FNO
FNO DSE
GINO
UPT
Heat-L-Sines
0.04
0.05
0.51
0.15
0.53
0.19
20.2
Wave-C-Sines
5.35
6.25
17.3
13.1
5.52
5.82
12.7
NS-Gauss
2.29
3.80
76.4
41.1
38.4
13.1
92.5
NS-PwC
1.58
2.03
33.1
26.0
56.7
5.85
100
NS-SL
1.28
1.91
18.8
24.3
29.6
4.48
51.5
NS-SVS
0.56
0.73
7.53
9.75
26.0
1.19
4.20
CE-Gauss
6.90
7.44
76.4
42.1
30.8
25.1
64.2
CE-RP
3,98
4.92
13.9
18.4
27.7
12.3
26.8
ACE
0.01
0.01
0.09
1.09
1.29
3.33
100
Wave-Layer
6.77
9.01
84.8
11.1
28.3
19.2
19.6
Poisson-Gauss
2.26
2.52
30.9
8.16
2.27
7.57
48.4
AF
1.00
1.09
10.1
4.48
1.99
2.00
45.7
Elasticity
4.31
4.63
11.9
5.53
4.81
4.38
12.6
Table 2: Benchmarks on datasets with Cartesian-grid space discretizations. Lowest and second lowest emrors highlighted. The errors for time-dependent problems are the lowest among three autoregressive strategies (see SM Section C.2p. All training datasets have 1024 trajectorics/samples.
Dataset
Median relative L' error [%]
(Cartesian grid)
RIGNO-18
RIGNO-12
CNO
scOT
FNO
NS-Gauss
2.74
3.78
10.9
2.92
14.2
NS-PwC
1.12
1.82
5.03
7.11
11.2
NS-SL
1.13
1.82
2.12
2.49
2.08
NS-SVS
0.56
0.75
0.70
0.99
6.21
CE-Gauss
5,47
7.56
22.0
9.44
28.7
CE-RP
3.49
4.43
18.4
9.74
31.2
ACE
0.01
0.01
0.28
0.21
0.60
Wave-Layer
6.75
8.97
8.28
13.4
28.0
Poisson-Gauss
1.80
2.44
2.04
0.68
11.5
grids. The resulting test errors in Table 2 show that RIGNO is very accurate Cartesian grids too and RIGNO-18 is the best performing model on out 9 datasets, often with several times lower errors than the baselines. We would also like emphasize that RIGNO is the only framework that displays acceptably low errors for every single dataset that we consider here, demonstrating its unprecedented robustness.
RIGNO scales with data and model size. An essential requirement for modern ML models is their ability to scale with data, i.e., the test error should reduce at a consistent rate as the number of training samples increases [14 training samples in Figure 2a observing that RIGNO-12 shows consistent scaling for the datasets. Another critical aspect of ML models their ability scale with model size 14 [3] Compared to data scaling, scaling with model size is more subtle as number of byperparameters influence model size and changing each of them, while keeping the others fixed, might lead to different outcomes. In SM Section F.3 we present detailed discussion and results on how RIGNO scales with respect to each of the hyperparameters that influences model size. Here, Figure 26 we only consider two (correlated) indications of model size, i.e., total number of model parameters and inference time. see from Figure 25 that RIGNO-12 consistently scales with respect to both the total number of parameters and inference run time.
To this aspect, we plot the error with respect the number of

---

Inference time [ms]
10
NS Gauss
2.68
4.12
13.5
20
NS-ST.
Error [%]
20
NS SVS
5
CE RP
Error [%]]
10
Prisemm-huaee
2
0.5
Lager
2-
256
512
1024
1024 Training dataset size
256
512
Training dataset size
Eluaticity
0.2
0.5
1.9 3.0
Parameters x10"
(a) Training dataset size.
(b) Model size (NS-PwC).
Figure 2: Scaling behavior of RIGNO with respect dataset and model size. All values correspond to autoregressive (AR-2, see SM Section C.2 test errors at ( 144
Resolution invariance and RIGNO. As mentioned before, a key requirement for a robust neural operator is to enforce appropriate form of resolution invariance 22. To test the abilities of RIGNO in this regard, we consider the NS-PwC dataset on random point clouds. The training resolution for RIGNO is 422. At inference, we test on this and several other resolutions, presenting the test errors in Figure 3a We see from this figure that the super-resolution test errоrs with the base version of RIGNO are very small and indicate that the proposed local aggregation scheme enables RIGNO to generalize well across higher spatial resolutions. The results of the same test on MeshGraphNet (see SM Section F.4) provide further evidence that local aggregations are essential for obtaining super-resolution invariance. Furthermore, the proposed edge masking technique allows us to enforce resolution invariance to a greater extent, particularly in the sub-resolution regime. Moreover, we observe that edge masking substantially improves the optimization process by allowing the model to learn weaker spatial correlations, while simultaneously cutting the computational cost by a factor of 2 when 50% of the edges are masked.
EM
L. 3
Ervor (96)
20
50% EM
Eror [96]
M
Error (X]
0.
16
32
42a 52" 64"76"
128
10 12
18 20
Inference resolution
Time index
Time index
(a) Edge masking (NS-PwC).
(b) Fractional pairing (NS-PwC).
(c) Time marching (NS-SL).
Figure 3: Ablation of the proposed training strategies. a) The effect of edge masking on resolution invariance. EM stands for edge masking. The y-axis shows the relative test error at = f14 when the spatial resolution of the inputs and outputs varies from the training resolution 642. b) Autoregressive test errors with and without fractional pairing fine-tuning. The initial solution (model input) is updated on vertical lines according to the AR-4 scheme (see SM Section C.2 These models are
exceptionally trained on snapshots with time resolution 4åt and up to time t16- c) Autoregressive
(AR-2) test errors with different time marching strategies. Time snapshots after the vertical line has not been seen during training and are considered as extrapolation in time.
On temporal processing with RIGNO. RIGNO can learn the solution operators of time-dependent PDEs and contains several elements in its design to enforce time continuity and enable accurate temporal approximation. To test these aspects, we focus on the NS-PwC and the NS-SL datasets. The underlying trajectories are available from the PDEgym collection [33 on a time interval of [o, 1]. sampled at equidistant times f, = nAt, with At 0.05 and 0 < < 20. We trained RIGNO with data sampled at times (ta, f4, fs. f12, )1a). The smallest time difference present in the training dataset is thus r = 4At. During inference, we unroll RIGNO using time steps (At, 2At, 34t, 4 according to the AR-4 scheme described in SM Section C.2 We compare the test errors at all possible time indices 1 < n < 20. Hence, we test the generalization capability of RIGNO to higher time resolution. The resulting test errors are shown in Figure (3b] With regular training, the errors with unseen time steps < 4At are significantly higher than with 4At, leading to the
8

---

observed jagged pattern in the plot of the error over time. The proposed fractional pairing fine-tuning approach effectively enhances the accuracy with small unseen time steps, leading to monotone error a
accumulation with & 2 14, where the output time is in-distribution.
We also consider time marching strategies arising from the choice of parameters o, B in (3). Here, we train and infer the model using the same time resolution (2At).. We present the results with all the three choices of time marching strategies, namely output, residual and derivative marching for the NS-SL dataset in Figure 3c We observe a slow error growth and very good extrapolation beyond the times seen during training (t > t14). Moreover, we see that the derivative time marching strategy with a = 1, 5 our evaluations.
is clearly superior in this example to the other two, leading us to prefer it during
Quantifying model uncertainties. Not only does edge masking enable RIGNO to be (approxi- mately) resolution invariant, but it also provides a simple method to quantify model uncertainties. As edge masking is random, invoking it during inference enables the model to provide multiple outputs for the same input. The mean of the outputs can be used as the model prediction as shown in SM Figure F.4 The same figure shows how the standard deviation of this ensemble of outputs is highly correlated with the test error. Thus, it can be used to provide an adaptive time marching strategy to reduce errors at inference, see SM Section for for details. In SM Section F. masking makes RIGNO more robust to noisy inputs as well.
we show that edge
4
Discussion
In this paper, we present RIGNO, an end-to-end GNN- neural operator for learning the solution operators of PDEs on arbitrary domains. Our model follows an encode-process-dee process-deodde paradigm [4]] Given a point cloud on which the input data is defined (the physical mesh), the encoder is a GNN that maps the input data into a latent space defined on a regional mesh, which can be interpreted as a downsampled version of the underlying point cloud. The input point cloud can be very general, as relational structure or graph must be specified on it. The processor maps the latent variables on the regional mesh into learned representations through multi-scale processing, implemented through multiple stacked pesssae-passinng neural networks. Finally, the decoder transforms the processed variables from the regional mesh to physical variables on the output point cloud. Again, the only graph in this layer connects the regional mesh nodes to the points in the output cloud. No relational structure or graph needs to be defined on the output point cloud, making RIGNO very general in terms of its abilities to handle data on arbitrary point clouds.
In addition to this basic construction (see Figure D. RIGNO introduces several novel elements to satisfy two important requirements. First, as desired in [16] 2 (approximate) resolution invariance is essential for structure-preseve in the regional mesh connecting points in the physical mesh and more importantly, edge masking, to enable approximate resolution invariance, feature which is hard to achieve for neural operators on arbitrary domains. Second, RIGNO is tailored process time-dependent PDEs. In this context, temporal continuity and control on eгror growth over time are key properties that the neural operator needs to uphold. Here, we have incorporated many elements in RIGNO to achieve these requirements, namely lead time conditioning, all2all training, and novel fractional pairing fine-tuning strategy. This allows us to maximize training steps for given number of trajectories. while allowing for general rollout strategies with arbitrarily small time steps at inference.
operator learning. RIGNO uses overlaps between the sub-regions
We have tested RIGNO on suite of challenging problems including PDEs with data struc- tured/unttructured meshes, random point clouds, and Cartesian grids. RIGNO is shown to be significantly more accurate than baselines designed for such arbitrary domains. We also observe that RIGNO is more accurate than oftte-of-theaan neural operators such as CNO, FNO and scOT, even for PDEs on Cartesian grids. The gains in accuracy are not marginal, and can be of an order of magnitude on many problems, demonstrating the expressive power of the proposed architecture. We attribute RIGNO's of-the-of accuracy and robustness primarily to the high expressive power of local message passing, which enables fineggraine manipulation and propagation of information. This advantage is especially prominent time-dependent and advection-doninate problems. For elliptic problems, where global interactions are more crucial, the gains are less pronounced. To mitigate the limitations of purely local interactions, RIGNO incorporates multi-scale edges in the processor, introducing long- range interactions that significantly improve accuracy without incurring notable
9

---

additional computational cost. Furthermore, we show that RIGNO scales with both model and data size. while having excellent resolution invariance and temporal continuity properties. It generalizes very well to interpolated and extrapolated time levels, and maintains its accuracy with very small time steps.
Related work. The popular MeshGraphNet [3 served as foundation for the formulation of RIGNO and thus admits basic architectural similarities. This approach constructs a graph on the input coordinates based on mesh and performs message passing steps directly on the physical coordinates. Yet, as confirmed by our empirical results (see SM Section F.4), this approach leads to strong dependence on the space resolution of the training data and cannot generalize to different resolution. Furthermore, the computational cost quickly becomes prohibitive as the number of mesh (input) points increases, hindering MeshGraphNet from scaling to larger problems. This issue is addressed in [10] by downsampling the initial mesh to a coarser mesh and performing more message passing steps on the coarse mesh. However, this approach still relies on the original mesh for resolving small scales and consequently inherits the dependence on the training resolution from MeshGraphNet. As evidenced by our experimental results, the multi-scale regional mesh construction of RIGNO mitigates the reliance of such graph-based approaches on the distribution of the input points and introduces the ability to scale to larger problem sizes (see Table F.2). Radius-based aggregation for projection onto a coarse grid has been explored in prior work and shown to be effective for non-graph-basee architectures. GINO [23] performs local aggregation of coordinates within fixed radius to map between an unstructured point cloud and a uniform mesh, enabling the use of FNO as the core processor. UPT employs transformers as the core processor, and thus does not require the mapped coordinates to lie on uniform grid. In this case, local aggregation is optionally applied in the encoder to compress the information from dense point clouds to a coarser representation. Both these methods rely on global fixed radius which is a hyperparameter and must be calibrated for each dataset. practice, the given point clouds are often more refined in regions with sharper gradients, making a fixed-radius approach destined to produce To address this, additional mechanisms such as limiting the number of incoming edges per node are necessary. By contrast, the adaptive strategy used in the current work determines aggregation radii based on local resolution, thereby mitigating both these issues while automatically ensuring full domain coverage without the need for manual calibration.
uneven number of edges per receiver node.
Limitations. This work focuses on PDEs two spatial dimensions. However, many of the considered datasets are comparable in complexity, grid resolution, and number of input-output pairs to typical 3D problems. In principle, RIGNO is readily extendable to PDEs in three dimensions, though the current paucity of publicly available 3D benchmarks limits our ability to explore this important direction. Furthermore, our focus here was on presenting the model and extensive empirical investigation. Important theoretical questions such as those of universal approximation of PDE solution operators with RIGNO will be considered in future work. Similarly, we have empirically demonstrated that RIGNO scales with the number of training samples; however, proving rigorous generalization error bound, such as that of [18], merits further analysis. Moreover, the scaling curves of Figure 2 suggest that large number of training samples are necessary for neural operators to achieve small error. This issue is common to all standalone neural operators
[33 and can be alleviated
Our results presented in SM Section F.9 show that RIGNO can largely benefit from transfer learning. Using RIGNO as the backbone for foundation model is another avenue for further exploration. Finally, downstream applications of RIGNO to tasks such as uncertainty quantification (UQ), inverse problems and PDE-constrained optimization will be considered in future work.
by building foundation models such those proposed in [K]
44 30)
Acknowledgements
The contribution of Siddhartha Mishra to this work was supported in part by the DOE SEA-CROGS project (DE-SC-002311)).
10

---

References
[1] B. ALKIN, A. FÜRST, L. SCHMID, L. GRUBER, M. HOLZLEITNER, AND J. BRANDSTET- TER, Universal physics transformrs: A framework for efficiently scaling neural operators, in The Thihty-eighth Annual Conference on Neural Information Processing Systems, 2024.
[2] F. BARTOLUCCI, E. DE BEZENAC, B. RAONIC, R. MOLINARO, S. MISHRA, AND R. ALAI- FARI, Representation equivalent neural operators: a framework for alias-free operator learning. Advances in Neural Information Processing Systems, 36 (2024).
[3] J. BRADBURY, R. FROSTIG, P. HAWKINS, M. J. JOHNSON, C. LEARY, D. MACLAURIN,
G. NECULA, A. PASZKF, VANDERPLAS, S. WANDERMAN-MILNN AND Q. ZHANG, JAX: composable transformations of Python+NumPy programs, 2018.
[4] J. BRANDSTETTER, D. E. WORRALL, AND M. WELLING, Message passing neural PDE solvers, in International Conference on Learning Representations, 2022.
[5] M. M. BRONSTEIN, J. BRUNA, T. COHEN, AND P. VELIČKOVIC, Geometric deep learning: Grids, groups. graphs geodesics, and gauges, arXiv preprint 1rXiv:210..1347, (2021).
[6] S. CAO, Choose a transformer: Fourier or galerkin, in 35th conference on neural information processing systems, 2021.
[7] T. CHEN AND H. CHEN, Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems, IEEE Transactions on Neural Networks, 6 (1995), pp. 911-917.
[8] L. EQUER, T. K. RUSCH, AND S. MISHRA, Multi-scale message passing neural PDE solvers, arXiv preprint arXiv:2302.03580, (2023).
[9] L. C. EVANS, Partial differential equations, vol. 19, American Mathematical Society, 2022.
[10] M. FORTUNATO, T. PFAFF, P. WIRNSBERGER, A. PRITZEL, AND P. BATTAGLIA, Multiscale meshgraphnets, arXiv preprint arXiv:2210.00612, (2022).
[11] Z. HAO, C. SU. S. LIU, J. BERNER, C. YING, H. Su, A. ANANDKUMAR, J. SONG, AND ZuIU, DPOT: auto-regressive denoising operator transformer for large-scale PDE pre- training. in Proceedings of the 41st International Conference on Machine Learning, ICML'24, JMLR.org, 2024.
[12] Z. HAO, Z. WANG, H. SU, C. YING, Y. DONG, S. LIU, Z. CHENG, J. SONG, AND ZHU, Gnot: A general neural operator transformer for operator learning, 2023.
[13] M. HERDE, B. RAONIĆ, T. ROHNER, R. KÄPPELI, R. MOLINARO, E. DE BÉZENAC, AND
S. MISHRA, Poseidon: Efficient foundation models for PDEs, arXiv preprint arXiv:2405.19101, (2024).
[14] KAPLAN, S. MCCANDLISH, HENIGHAN, B. BROWN, B. CHESS, R. CHILD, S. GRAY. A. RADFORD. J. Wu, AND D. AMODEL, Scaling Laws for Neural Language Models, Jan. 2020. arXiv:2001.082 [es, stat].
[15] N. KOVACHKI, S. LANTHALER, AND S. MISHRA, On universal approximation and enror bounds for fourier neural operators, Journal of Machine Learning Research, 22 (2021), pp. Art- No.
[16] N. KOVACHKI, Z. LI, B. LIU, K. AZIZZADENSHELI, K. BHATTACHARYA, A. STUART, AND A. ANANDKUMAR, Neural operator: Learning maps between function spaces, arXiv preprint arXiv:2108.0844 1v3, (2021).
[17] R. LAM, A. SANCHEZ-GOnE GONZALEZ, M. WILLSON, P. WIRNSBERGER, M. FORTUNATO,
F. ALET, S. RAVURI, EWALDS, Z. EATON-ROSEN, W. HU, A. MEROSE, S. HOYER, G. HOLLAND, O. VINYALS, J. STOTT, A. PRITZEL, S. MOHAMED, AND P. BATTAGLIA, Learning skillful medium-romee range global weather forecasting, Science, 382 (2023), pp. 1416 1421.
11

---

[18] S. LANTHALER, S. MISHRA, AND G. E. KARNIADAKIS, Error estimates for DeepONets: A deep learning frameworr in infinite dimensions, Transactions of Mathematics and Its Applica- tions, 6 (2022), p. tnac001.
[19] Z. LI, D. Z. HUANG, B. LIU, AND A. ANANDKUMAR. Fourier neural operator with learned deformations for pdes general geometries, Journal of Machine Learning Research, 24 (2023), pp. 1-26.
[20] Z. LI, N. B. KOVACHKI, K. AZIZZADENESHELL, B. LIU, K. BHATTACHARYA, A. STUART, AND A. ANANDKUMAR, Fourier neural operator for parametric partial differential equations, in International Conference on Learning Representations, 2021.
[21] Z. LI, N. B. KOVACиКI, K. AZIZZADENESHELI, B. LIU, K. BПATTACHARYA,, A. M. STU- ART, AND A. ANANDKUMAR, Neural operator: Graph kernel network for partial differential equations, CoRR, abs/2003.0344 (2020).
[22] Z. LI, N. B. КОVАСиıКı, K. AZIZZADENESHELLI B. LIU, A. M. STUART, K. BIAT-- TACHARYA, AND A. ANANDKUMAR, Multipole graph neural operator for parametric partial differential equations, in Advances in Neural Information Processing Systems (NeurIPS), H. Larochelle. M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, eds.. vol. 33, Curran Asso- ciates, Inc., 2020, pp. 7766-
[23] Z. LI, N. B. KOVACHKI, C. CHOY, B. LI, J. KOSSAIFI, P. OTTA, M. A. NABIAN, M. STADLER, C. HUNDT, K. AZIZZADENESHELI, AND A. ANANDKUMAR, Geometry- informed neural operator for large-scale 3d PDEs, in Thirty-seventh Conference on Neural Information Processing Systems, 2023.
[24] L. E. LINGSCH, M. Y. MICHELIS, E. DE BEZENAC, S. M. PERERA, R. K. KATZSCHMANN, AND S. MISIRA, Beyond regular grids: Fourier-based neural operators on arbitrary domains, in Proceedings of the st International Conference on Machine Learning, R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp, eds., vol. 235 of Proceedings of Machine Learning Research, PMLR, 21-27 Jul 2024, pp. 30610-30629.
[25] Z. LIU, H. HU, Y. LIN, Z. YAO, Z. XIE, Y. WEI, J. NING, Y. CAO, Z. ZHANG, L. DONG, ET AL.. Swin transformer V2: Scaling up capacity and resolution, in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2022, pp. 12019.
[26] Z. LIU, H. MAO, C.-Y. WU, C. FRICHTENHOFER, T. DARRELL, AND S. XIE, A ConvNet for the 2020s, in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2022, pp. 11976-11986.
[27] I. LOSHCHILOV AND F. HUTTER, Decoupled weight decay regularization, in International Conference on Learning Representations, 2019.
[28) L. LU, P. JIN, G. PANG, Z. ZHANG, AND G. E. KARNIADAKIS, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence, 3 (2021), pp. 218-229.
[29] A. J. MAJDA, A. L. BERTOZZI, AND A. OGAWA, Vorticity and incompressible flow, Appl Mech. Rev., 55 (2002), pp. B77-B78.
[30] M. MCCABE, B. R.-S. BLANCARD, L. H. PARKER, R. OIIANA, M. CRANMER, A. BIETTI,
M. EICKENBERG, S. GOLKAR, G. KRAWEZIK, F. LANUSSE, M. PETTEE, T. TESILEANU,
K. CHO, AND S. Ho, Multiple physics pretraining for spatiotemporal surrogate models, in The Tiihth- Annual Conference on Neural Information Processing Systems, 2024.
[31] S. MISHRA AND A. E. TOWNSEND, Numerical Analysis meets Machine Learning. Handbook of Numerical Analysis, Springer, 2024.
[32] S. MOUSAVI, S. WEN, I.. LINGSCH, M. HERDE, B. RAONIĆ, AND MISHRA, Datasets for RIGNO, 2025.
12

---

[33] T. PFAFF, M. FORTUNATO, A. JANCHEZ-GOZZANN AND P. BATTAGLIA, Learning mesh- based simulation with graph networks, in International Conference on Learning Representations. 2021.
[34] M. PRASTIIOFER, T. DE RYCK, AND S. MISIIRA, Variable input deep operator networks, arXiv preprint arXiv:2205. 11404, (2022).
[35] A. QUARTERONI, A. MANZONI, AND F. NEGRI, Reduced basis methods for partial differential equations: introduction, vol. 92, Springer, 2015.
[36] A. QUARTERONI AND A. VALLI, Numerical approximation of Partial differential equations. vol. 23, Springer, 1994.
[37] P. RAMACHANDRAN, B. ZOPH. AND Q. V. LE, Searching for activation functions, arXiv preprint arXiv:1710.05941, (2017).
[38] B. RAONIĆ, R. MOLINARO, T. DE RYCK, ROHNER, BARTOLUCCI, R. ALAIFARI,
S. MISHRA, AND E. DE BÉZENAC, Convolutional Neural Operators for robust and accurate learning of PDEs, Dec. 2023. arXiv:2302001178 [cs].
[39] A. SANCHEZ-GONZALER J. GODWIN, T. PFAFF, R. YING, LESKOVEC, AND
P. BATTAGLIA, Learning to simulate complex physics with graph networks, in International conference on machine learning, PMLR, 2020, pp. 8459-8468.
[40] A. SANCHEZ-GONZALEZ N. HEESS, J. T. SPRINGENBERG, J. MEREL, M. RIEDMILLER, R. HADSELL, AND BATTAGLIA, Graph networks as learnable physics engines for inference and control, in International Conference on Machine Learning, PMLR, 2018, pp. 4479.
[41] L. SERRANO, L. LE BOUDEC, A. KASSAI KOUPAI, T. X. WANG, Y. AND P. GALLINARI, Operator learning with neural fields: Tackling pdes on general geometries, in Advances in Neural Information Processing Systems, 2023,
VITTAUT,
[42] L. SERRANO, T. X. WANG, E. LE NAOUR, J.-N. VITTAUT, AND P. GALLINARI, Aroma: Preserving spatial structure for latent pde modeling with local neural fiels, in Advances in Neural Information Processing Systems, 2024.
[43] N. SRIVASTAVA, G. HINTON, A. KRIZHEVSKY, SUTSKEVER, AND R. SALAKHUTDINOV, Dropout: simple way to prevent neural networks from overfitting, The journal of machine learning research, 15 (2014), pp. 1929-1958.
[44] S. SUBRAMANIAN, HARRINGTON, K. KEUTZER, W. BHIMJI, D. MOROZOV, M. W. MAHONEY, AND A. GHOLAMI, Towards foundation models for scientific machine learning: Characterizing scaling and transfer behavior, in Advances in Neural Information Processing Systems, A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, eds., vol. 36,
Curran Associates, Inc., 2023, pp. 71242-71262.
[45] H. WU, H. LUO, H. WANG, J. WANG, AND M. LONG, Transolver: A fast transformer solver for pdes general geometries, in International Conference on Machine Learning, 2024.
13

---

Supplementary Material
RIGNO: A Graph-based Framework for Robust and Accurate Operator Learning for PDEs on Arbitrary Domains
A Graph Construction
A.1 Encoder and decoder edges
The edges between the physical and regional nodes in the encoder and in the decoder are defined based on support neighborhoods around each regional node. We consider circular neighborhood around each regional node (see Figures [ and la) and calculate the radii of the circles such that the union of all the circles covers the whole domain. To achieve this, we first construct a Delaunay triangulation on the regional nodes. We then calculate all medians of all the simplices in the triangulation. The minimum radius of node then calculated as two-thirds of the longest median that crosses that node. The minimum radius is then multiplied by overlap factor, which is hyperparameter of the model. A similar method can be applied to three-dimensional point clouds as well. Each regional node is then connected to all the physical nodes that lie inside its corresponding support sub-region. With a larger overlap factor, the regional nodes receive information from more physical nodes in the encoder. Too large overlap factors might be unnecessary since regional nodes are assumed to encode information locally. However, in the decoder, overlap factor plays a more important role. With minimum overlaps, many physical nodes will end up being connected to only single regional node, which might result in overly biased decoded latent features. Therefore, we generally use different overlap factors in the encoder and in the decoder.
A.2 Regional nodes
By defining the edge connections between the regional and the physical nodes based on support sub-regions in the encoder and in the decoder, RIGNO allows for unstructured physical and regional meshes. Within this framework, even for an unstructured set of physical nodes, the regional nodes can still follow structure (e.g., lie completely arbitrary. However, constructing an unstructured regional mesh is favorable since local resolutions of the physical mesh can be inherited by the regional mesh. To this end, we always define the regional nodes by randomly sub-sampling the given physical nodes by a factor that is given as hyperparameter. The set of regional nodes may be different during inference or in each training epoch.
uniform grid) since the choice of the regional nodes is
The sub- sub-sampling factor is fixed hyperparameter, which makes the resolution of the regional mesh dependent on the resolution of the physical mesh. If model only trained with a fixed mesh resolution, changing the resolution at inference might affect its accuracy. However, since the choice of the regional mesh nodes is completely arbitrary, the resolution of the regional mesh can be corrected in order to correspond to the resolution of the training regional meshes, For zero-shot super-resolution inference, this can be easily done by sub-sampling the physical mesh with a higher sub-sampling factor. For zero-shot sub-resolution inference, however, a suitable mesh-refinement algorithm must be employed. In this work, we refine the mesh by first constructing a standard triangulation, randomly picking the necessary number of triangles, and adding the centroid of those triangles to the regional mesh. With this method, we ensure that local resolutions are inherited from the original mesh.
A.3 Processor edges
The processor of RIGNO consists of sequence of message-passinn blocks with skip connections, which act on the regional nodes and the edges that connect them. We define the first level of the bidirectional edges between the regional nodes as all the unique pairs of connected vertices in Delaunay triangulation on the regional nodes. In addition to the first level edges, we also define edges that connect long-distan distance regional nodes, which allows learning multiple spatial scales and enables a more efficient propagation of information in a message-passing step. The construction of these edges is illustrated in Figure [] as multi-scale edges. For the second level, we randomly choose subset the regional nodes and repeat the above process on them. We repeat this process multiple times by sub- samplinm the subsets again.
14

---

The structure of the regional mesh can be regarded as a hierarchy of multiple sub-meshes with different resolutions and shared nodes (see Figures Figures/A.I and the multi-scale edges box of Figure In one message-pasae passing step, the coarser sub-meshes, which have longer edges, allow the information to flow in larger scales, while the finer sub-meshes allow for assimilating small-scale details. All regional nodes are part of the finest sub-mesh, while only few nodes take part in the coarsest one. a single message-passing step, information from longer distances. Because of this, is essential to have multiple message-passing steps in the processor to allow for the information to propagate through the whole domain.
not possible for all the nodes to receive
A.4 Fixed and variable physical meshes
For some applications engineering (e.g., shape optimization), the physical mesh changes for every sample in the training dataset or at inference. Since the coordinates of the input discretization will change, the graph construction needs to be carried out every time. For training with large datasets, these graphs can be constructed only once as pre-processing step. At inference, however, the time for constructing the graph is considerable compared to the execution of the model. On the other hand, another class of applications in science and engineering only require training and inference on fix unstructured physical mesh. For these problems, constructing the graph can only done once.
For the latter class of problems, RIGNO is trained using fixed physical mesh. However, if the graph construction is only done once during training, regional mesh will consequently be fixed too. this case, although the encoder and the decoder are invariant to the input discretization, the processor might be biased to that specific fixed regional mesh. Such a model can exhibit poor performance when it is inferred with a different regional mesh. In order to avoid this situation, we update the regional mesh before every epoch during training.
(a) Aggregation area of regional nodes.
(b) Level 1 regional mesh edges.
(c) Level 2 regional mesh edges.
(d) Wrapped regional mesh edges.
Figure A.1: Construction of regional mesh edges and support sub-regions with periodic boundary conditions. The dashed lines show the domain boundaries. Large and small purple dots represent the main regional nodes and the ghost regional nodes, respectively. Black lines represent the main edges, and the purple lines the cross-boundary edges. Note that in Figure Ia the parts of the domain that are close to the domain boundaries and do not lie in any circle are covered by at least one regional node on the other side of the boundary.
A.5 Periodic boundary conditions
The majority of the problems considered in this work have periodic boundary conditions. The flexibility graphs allows us to enforce periodic boundary conditions in the architecture without any significant computational overhead. To achieve this, we make the following three changes to our framework for problems with periodic boundary conditions: 1) Initialize the node features with a periodic encoding of the coordinates instead of raw coordinates; 2) extend the domain before defining the encoder and the decoder edges; and 3) extend the domain before defining the processor edges.
In order to motivate the above changes, let us consider two nodes, one close to the left and one close the right boundary. With (normalized) raw coordinates, the initial features of these two nodes will be very different as one will take value close to and the other one takes value close to 1. With periodic boundary conditions, however, these two points will be close to each other if frame of the reference domain is slightly shifted. The invariance to shift in the frame of reference means that these two nodes should be connected in the graphs. simpler words, we assume during graph construction that the both ends of domain with periodic boundary conditions connected. In
15

---

one-dimensional settings, this set-up corresponds to considering circular domain, see Figure A.2 for an illustration. For square two-dimensional domain, connecting the two ends in both direction corresponds to torus.
Inspired by the ghost node technique in classical numerical methods, we extend the domain in all directions by constructing a x tiled structure, where each tile contains the repeated coordinates of the given point cloud. All the Delaunay triangulations in the graph construction phase are built based on the extended domain. Figure A. shows the resulting support sub-regions, as well as the extra edges that are introduced because the domain extensions for a minimal set of regional nodes. We refer to these extra edges as cross- moos--bundary edges. Notice how all the white areas in Figure A.la are covered by at least on circle on the other side of the domain. In order avoid introducing extra costs to the architecture, the ghost nodes are discarded after graph construction and only the cross-boundary edges are kept. order to keep the extra edges, we replace each ghost node with its corresponding node in the original point cloud. The result is illustrated in Figure A.
Plynioal regional olgas (Encader)
Regional Mesh edges (Processor)
x 1,
Regional plastcaa odges (Decsder)
Edge morkin with prsbubility 0.5
Figure A.2: Graph structure for periodic boundary conditions in one-dimensional settings. The green dots represent physical nodes, the purple circles represent regional nodes, and the solid lines represent edge connections. Note that only a few physical nodes are visualized. The bottom right part of the figure partly shows the receptive fields of two regional nodes in the encoder. Here. the red lines show the edges that are randomly masked, and hence will not be used in that particular message passing block.
A.6 Structural features
The corresponding coordinates of the nodes are pre-processed to provide the initial structural features of the graph entities. In order to impose translation equivariance to the model, we consider linear normalized coordinates Ç: E [-1, +1] and angular normalized coordinates a, E [0, 2T] defined as
S: : 2(x;- x,)LL 1 and a; T(. + 1). The structural features of edge eij connecting node i node j sp.j := Kč. [12): with S.j == (C, (2V/))))) being the relative coordinates of the connected nodes normalized by the diameter of the domain, 2/d. We denote the vertical concatenation of column vectors by o. The structural features of node v, with general boundary conditions are simply its normalized coordinates P := Si With periodic boundary conditions (see Section A.5 we concatenate sines and cosines of angular normalized coordinates with Ktreq frequencies: 9 [sin(ka;) cos(ka;)) for 1 < k < Kireq-
16

---

B Details of the GNN Components
We consider two underlying sets of nodes (see Figure V2 = {v denotes the physical nodes, which corresponds to the unstructured point cloud which the input (and output) data is defined. The regional nodes are denoted by 9 = {0, Ji and the bidirectional edges that connect them by
Ex {e (e Jay- Ex = {e 2-+ hij and Ex- = {ee) + his denote the unidirectional
edges the encoder and the deocder, respectively.
The first step is to embed the initial features all nodes and all edges (see Section A.6) to suitable high-dimensional features. The other model inputs (initial condition, spatial parameters, and time) only embedded the features of the physical nodes. For regional nodes, we also embed the radius of their support sub-region Rg. This done by passing all the initial features through the embedder feed-forward blocks (FF, see Figure
(.(0)
FF
embedder (u, Pi t,,
(Physical nodes),
0
FFR
embedder (4,: R&)
(Regional nodes),
embodder (4i;))
(Physical to regional edges),
(5)
(t
embodder (pts))
ot
(Regional mesh edges),
St (t))
2- embodder (4:ss)
(Regional to physical edges).
Next, the encoder applies single message passing step to the embedded features of the physical nodes and transmits the information to the regional mesh,
(-+..(()
FT-0 - P.(O)
(0)
2-9H.(0)
m ij
encoder (v)
eij
(1)
(R)(0)
F
(0)
Σ
(--R.(0)
v;
encoder
M
mji
(6)
JEN,
v U.
(1)
3 (.)
FT encoder
B00))
The processor passes this latent information on the regional mesh through P sequential message- passing blocks with residual updates and skip connections. The p-th message-passing step in the processor updates the regional nodes and the regional mesh edges
Ap m 1j (x+11) U
(72 (u: U,
H.(p) (..))
(o-1), 1)),
eij
=0,
F1 st(e)
processor (V
W.
m
(p-)),((-11),
(7)
JEM.
(p)
(x-1)
(x-1)
ij
eij
m ij
The decoder then transfers processed information from the regional mesh back to the physical mesh by performing single message passing step,
n 90)3((0)
FF .P.
+1)
(,(1)
m
decoder (v,
. C,j
(2) (1)(1)
2F
(1)
m 90FF1(0)
(8)
v Vi
decoder
Mi
JEN.
Finally, in the second step of the decoder, the feature vectors of the physical node are independently scaled down to the desired dimension to form the output
Fooutput (v) (2.(2))
(9)
17

---

C Temporal Processing
C.1 Time stepping strategies
As described in the main text, we consider three time stepping strategies, namely direct outpur stepping, residual stepping, and time-derivative stepping. It is important to normalize the outputs correspondingly in each of these strategies since the nature and the statistics of the target variables are different in each case. Therefore, we pre-compute the appropriate statistics in each strategy based on the training samples and use them normalize and de-normalize the inputs and the outputs. With direct output stepping, the statistics of the inputs and the outputs are trivial. For residual stepping. pre-compute the statistics of the residuals with all the lead times that are present in the training. For time-derivative stepping, we similarly consider all lead times and compute the statistics of the residual divided by the lead time. Figure C.1 illustrates flowcharts for input and output normalization with each stepping strategy.
ronidaal
dertvatroe
EOss
rotraated
rosiral
edlraated roiiul
ralreated
entonainn
largst
roidual
(a) Output stepping.
(b) Residual stepping.
(c) Time- dime-derivative stepping.
Figure C.I: Flowcharts of input and output normalization strategies. Blue boxes represent raw data and green boxes represent normalized data. We apply a min-max normalization for the ipput time which is omitted in the flowcharts for simplicity. The loss is computed on normalized values while the errors are computed on the raw targets and de-normalized estimates. Note that the raw target and the de-normalized estimates may be normalized in different manner or with different statistics (e.g.. global equation- aquation-wide statisties as in 13 before the error is computed on them.
C.2 Temporal inference
Given any set of partial lead times [T. that sum up to a time E Qr time-continuous operator Go can be applied sequentially to reach the solution time t:
utt)= Go (Go (gp ....c,
n-3
n 2
Σ
To T
c,
Tss Fn-1
c,
Tea Tn
(10)
This gives an infinite number of ways to infer the learned operator for reaching u(t). However, since Go is not exact, the approximation errors generally accumulate and can lead to very high errors for long-time estimations. With larger lead times it is possible to decrease the number of times that the operator is applicd to completely avoid autoregressive inference by directly using = t.
In all the time-dependent datasets considered in this work, for simplicity of exposition, we have considered N₂ 20 or N, = 21 time snapshots {t., N. Unless otherwise stated, only snapshots with even indices up to f14 are used for training. Snapshots with > f14 are reserved for time extrapolation tests, and snapshots with < (14 with even indices are reserved for time interpolation tests. Error at snapshot (14 is used as a reference for the benchmarks. Since the trained operator is continuous in time and can be inferred with different lead times, there are infinitely many ways to reach the solution at f14 via autoregressive inference of the trained operator. The optimal autoregressive
18

---

path can be different for different operators. Figure C.2 illustrates the training and test time snapshots as well as different autoregressive paths for reaching 114- In the presented benchmarks in Tables
1 and 2 we calculate the error with 3 different autoregressive paths: i) (T: 2 which is denoted by AR-2; ii) (π = 4âL, 72 = 4A1, 73 = 4A1, 74 = 2 which is denoted by 4R-4; and iii) { 14 denoted by DR. We always report the lowest error among these three rollout paths.
Tix AA0
A
Used for training
Final in-distribution snapshot
Used for time-interpolation tests
Final extrapolation snapshot
Used for time-cxtrapolation teats
Figure C.2: Snapshots used for training, time interpolation, and time extrapolation. The in-distribution and extrapolation evaluation snapshots are marked with thick borders. The lower arrows show the autoregressive steps for evaluating a learned operator at t1 with Tmax = this rollout path, lead time (At) that is smaller than the time resolution of the training dataset (2At) is applied 14 times in the evaluation. The upper arrows show the autoregressive steps for evaluating a learned operator 220 with mmax = 4At (5 steps).
At (14 steps). Note that in
C.3 Fractional pairing
We use the fractional pairing strategy to achieve a truly coneinontinuuus operator that can be inferred with lead times smaller than the smallest time difference present the training dataset. In fractional pairing, we approximate the solution at an intermediary time by inferring a model that has been trained on the original time resolution. This approximation is then given as input to the same model. By choosing a suitable intermediary time, we can have an input-output pair that is not present in the training dataset and have a smaller lead time than the time resolution of the training dataset. This idea is illustrated in Figure C.3 With this, we can either calculate the loss gradients for both of the inferences, or we can cut the gradient at the intermediary time and treat the intermediary solution In the
as a normal model input, which is referred to as GradCui (see bottom left of Figure experiments of this work, we always use fractional pairing with GradCut.
If the approximation at the intermediary time is not accurate enough, fractional pairing strategy can lead to training instabilities. Therefore, it recommended use this technique only as a fine-tuning stage. Furthermore, since the approximations are made with fractional lead times, i.e., lead times that are not present in the training dataset, is also essential to employ other mechanisms that help generalize to fractional lead times with regular training strategies, such as consistent time marching strategies and conditiooe normalization layers. further consideration that helps the stability of this strategy is avoiding out-of-distribution lead times in the first forward pass, i.e., avoiding lead times that are smaller than the time resolution of the training dataset. Such small lead times can be seen as an extrapolation in the seen distribution of lead times during training. It has been observed that even without fractional pairing, the models generalize better to interpolated lead times than they do to extrapolated ones. Therefore, we can approximate the intermediary solution with an interpolated lead time, and use a smaller lead time only in the second forward pass. In Figure C this corresponds to only using the two right-most additional pairs in the green box. Note that none of the additional pairs of the orange box will be viable in this setting.
D Datasets
RIGNO is tested on suite of 13 datasets 32 which covers operators arising from compressible Fuler (CF) equations of gas dynamics, incompressible Navier-Stokes (NS) equations of fluid dynamies, wave equation (WE), heat equation (HE) of thermal dynamics, Allen-Cahn equation (ACE), Poisson's
19

---

GI GRADDENT FLOW
GrotCaa
aupplact triiting datact
Medel-infirmol srapulat
Pairaleainal by allbil maining
Pairobesinal by frsctioaa pairing
Figure C.3: Schematic of the fractional pairing strategy. The additional pairs obtained by fractional pairing correspond to smaller lead times than the smallest time difference in the original trajectory. On the bottom left of the figure, the forward passes (solid black arrows) and the backward gradient flows (dashed arrows) are illustrated for an existing pair (f,, f,). With the GradCut setting, the gradient flow is stopped at t..
equation, and force balance equation in solid mechanics (FB). Table D.I summarizes the main characteristics of all the datasets. A random sample of each dataset is visualized in Figures D.2 D.I4 This collection of datasets covers diverse multi-scale, nonlinear, complex phenomena systems arising from PDEs such shock waves, shear layers, mixing, and turbulence. This makes effort to learn the underlying solution operator from data worthwhile for these datasets; since running simulations of these systems with any classical numerical method would require extensive computational resources. Moreover, by considering these datasets, we cover diverse set domains and meshes from Cartesian uniform grids on unit square to structured and unstructured meshes fixed and varying irregular geometries.
The problems in the 8 depeedent datasets from 13 are all solved in the unit spatio-temporal domain Q2 x 22 [o. T] [o, 1]3. The final time is T except for the ACE dataset where T = 0.0002. For each dataset, a parameterized distribution of initial conditions considered, and initial conditions are randomly drawn by sampling those distributions. The distribution of the initial conditions characterizes the dataset, the level of complexity of the dynamics captured the phenomenological features exhibited the solution function, and the corresponding solution operator. In all these datasets, N, 21 snapshots of each sample are collected at times 2 {t,, n At == where N₂ = 20. Each snapshot originally has 128x 128 resolution, with uniform discretization in random point cloud with these datasets, randomly shuffle the coordinates and only keep the first 9216 coordinates in our experiments. Figure D.I illustaats a few random point clouds obtained from this procedure.
N,-1 with uniform spacing At; except for the ACE and the Wave-Layer datasets
both spatial directions. In order to test the baselines
Figure D.I: Random point clouds (9216 points) obtained from 128x uniform grid.
20

---

Table D.1: Abbreviations and main characteristics of datasets. IC: initial condition; BC: boundary
conditions.
Abbreviation
PDE
BC
Domain (varies)
Time
Source
Figure
Heat-L-Sines
HE
Dirichlet
L-shape (X)
D.2
Wave-C-Sines
WE
Dirichlet
circle (X)
D.3
NS-Gauss
INS
periodic
square (X)
NS-PwC
INS
periodic
square (X)
D.5
NS-SL
INS
periodic
square (X)
D.6
NS-SVS
INS
periodic
square (X)
D.
CE-Gauss
CE
periodic
square (X)
D.8
CE-RP
CE
periodic
square (X)
D.9
ACE
ACE
periodic
square (X)
Wave-Layer
WE
absorbing
square (X)
Poisson-Gauss
PE
Dirichlet
square (X)
X
13
D.13
AF
CE
various
airfoil (V)
X
Elasticity
FB
various
square with holes (V)
X
19
D.1 Equations
D.1.1
Incompressible Navier-Stokes equations (INS)
The incompressible Navier-Stokes equations are given by
V
(Conservation of mass),
(11)
0,0 + (v V v Vp+ (Conservation of momentum).
where v : Ω4 x 22p + I² is the velocity field with components (ea, vy) "y 92 x Q2 - R is the pressure, and is the fluid viscosity. In all the INS datasets, periodic boundary conditions considered and small viscosity v = 4 x 10-4 is applied to the high Fourier modes only.
D.1.2 Compressible Euler equations (CE)
The compressible Euler equations in the conservational form are given by
0,u +
S
u (p. pv, E) T)T,
F (pu,p Pn + (E p)
(12)
E p|lll1?
7
OΦ Os
(0,0, - -an--p-pp ppuy) O6
S (0, -p- 0, pv
Oy
with - = 1.4, where F is the flux matrix, E is the total energy density (energy per unit volume), and is the gravitational potential field. The gravity source term not considered in any the datasets of this work.
D.1.3 Allen-Cahn equation (ACE)
The Allen-Cahn equation describes phase transition in materials science and is given by dyu - 2 = 22 - -w),
(13)
where e is the reaction coefficient.
D.1.4
Wave equation (WE)
The wave equation is considered with fixed-in-time propagation speed that depends on the spatial location. This setting describes the propagation of acoustic waves in inhomogeneous medium:
du v2 0.
(14)
21

---

The wave equation can be written in the form of (T) by considering the first time derivative 0 as variable and writing as system of equations. The propagation speed can either be added as an input to the differential operator, or as another variable with 0,c = 0. In the former case, it should also be added to the inputs of any solution operator.
D.1.5
Heat equation (HE)
The heat equation describes the diffusion heat in medium and is given by
V²v²u,
(15)
where a the diffusivity constant.
D.1.6 Poisson's equation (PE)
Poisson's equation is given by
Vu f.
(16)
where f is given source function and u the solution.
D.1.7 Force balance equation (FB)
The force balance equation of solid body is given by
pd + V o=
(17)
where p denotes mass density, u the displacement field, and the stress tensor. The strain field is related to the displacement field via kinematic relations. The system is closed by considering constitutive model which describes the relation between strain and stress.
D.2
Dataset details
D.2.1
Heat-L-Sines
This dataset considers time-dependent solutions of the heat equation over an L-shaped domain which is formed by removing the rectangular portion [1, 1.5] x [0.5, I] from the larger rectangle [0.5, 1.5] x [0. 1]. We consider constant diffusivity = 1.0 and zero homogeneous Dirichlet boundary conditions. We prescribe superposition of sinusoidal modes the initial condition:
sin(Tmrii) sin (TT(22))))
uo 121, 22, 11)
Σ pm
(18)
0.5
where μ []1, M 10, and each Pm is sampled independently and uniformly from [-1, 1]. A Delaunay- triangular mesh with characteristic length 0.008 is used, yielding about 14,047 nodes and 27,590 elements. We simulate the heat equation for 100 time steps with a time step size of At = 5 x 10-5 and keep 21 (including the initial condition) equally spaced snapshots from each trajectory.
t4
ts
t= f12
t16
Figure D.2: Snapshots of solution trajectories of the Heat-L-Sines dataset.
22

---

D.2.2 Wave-C-Sines
This dataset pertains to time-dependent solutions of the two-dimensional wave equation over circular domain radius 1, centered at (0.5, 0.5). We set 2.0, constant time and space, and impose homogeneous Dirichlet boundary conditions (u 0) on the boundary of the circle. The initial condition for u at is specified as a superposition sinusoidal modes:
M
uo (2. v) =
Σ 1
Ois (1+ sin(rix) sin(Tjyy)
(19)
M
where K = 10, r 0.5, and aij ([111-- We employ Delaunay-based triangular mesh with characteristic length 0.015, resulting in roughly 16,431 nodes and 32, 441 elements. A time step of Al = 0.001 is used, and solutions are computed for 100 time steps via a finite element solver. We then downsample the computed solution fields in time, keeping 21 (including the initial condition) equally spaced time frames, which constitute the dataset.
to
t= t
t= Eg
t= t1
t= t16
1e-2
Figure D.3: Snapshots of solution trajectories of the Wave-C-Sines dataset.
D.2.3 NS-Gauss
The initial velocities are obtained from an initialized vorticity defined by w(x, y) = 0, U, u. by
using the incompressibility condition 29 of M 100 Gaussian distributions:
The vorticity field initialized linear superposition
(x
x2)² 2 + (y y1)
(x,
exp
(20)
σi
20
with 04 ~ 1.11) o; Up. 10.01.0.11,- 3t Mo.11. and y: 11. 11- A random sample of this dataset is
visualized Figure D.4
t= t4
t= te
t= t12
t= t16
Figure D.4: Snapshots of solution trajectories of the NS-Gauss dataset. Each row corresponds to one variable.
D.2.4 NS-PwC
The initial velocities are obtained from an initialized vorticity by using the incompressibility condition. The initial vorticity " of the flow field is initialized collection of piece-wite wise constant squares
w (x y) Cij in [1]1)22] X 13 Vy].
(21)
23

---

with 24 Bi for i = 0, 1,,,,..... P. and Ct ~ 1,11 The number of squares in each direction
P 10, resulting in 100 squares in total. random sample this dataset is visualized in Figure
D.5
t= to
t= t4
t= t8
t=
t12
t t16
-0.5
0.5
0.0
-0.5
Figure D.5: Snapshots of solution trajectories of the NS-PwC dataset. Each corresponds variable.
D.2.5 NS-SL
The velocity field v = (Uz, 02) is initialized as
tanh 2n 0-225
for y + 0z(2) S otherwise
(2r, v)
tanh 2n 0.75-v
(22)
0 y) O,
where as: [0, R perturbation of the initial data given by
(8(x)= = E a sin sin(2wkr - Bx))
(23)
The parameters sampled p 117.1..10..11.23- 0g Uj, Mj.22)) 5 0.025.
Mo..8.0.12]. and & ~ U 0.000.623.006255) A random sample of this dataset is visualized in Figure
D.6
t= to
t= t4
t= t8
t= t12
t= t16
Figure D.6: Snapshots of solution trajectories the NS-SL dataset. Each row corresponds to one variable.
D.2.6 NS-SVS
The vorticity of the flow is initialized as
w (z.v) 61(2, y) w w'(x (x,y).
(24)
24

---

where
w (x,y) 8(x r
dr,
(25)
T2
|
(x,
(26)
(r) = 80
(r + 1) - 4(r 2 2) + 6r - 4(r
2 2 + (r 1)
(27)
7m
10
1 + +0.2sin(2Tz)
r (( E T 1 y 2
+
sin (2 + B.)))
(28)
The random variables a; and 3 are given by a ~ 1. 0....0003125) 8 ~ 1)
The parameter p is fixed
as p
128
5
A random sample of this dataset is visualized in Figure D.7
t= to
t= t4
t= t8
t= t12
t= t16
0.25
0.00
-0.25
0.25
-025
Figure D.7: Snapshots of solution trajectories of the NS-SVS dataset. Each row corresponds to one variable.
D.2.7 CE-Gauss
The w field is initialized as linear combination of M 100 Gaussians
M
(r r; )2 + (v y1)
w (x,
exp
(29)
20
where a, ~ U(-1,11 ~ Mp.o1...1).. ~ Mo.1j. and ys ~ Up.1). The density and pressure are initialized with p == 0.1 and 2.5, both constant in the whole domain. Similar to the NS-Gauss dataset, the initial velocity field is obtained from the initialized vorticity by using the incompressibility condition. random sample of this dataset is visualized in Figure D.8
D.2.8 CE-RP
This dataset corresponds to the four-quadrant Riemann problem, which generalizes the standard Sod shock tube to two spatial dimensions. The domain is divided into a grid of P x P square subdomains Di,; {ax, y) T 1-1 < x < p. <y< < $1. The solution on each of these subdomains is constant and initialized as
j-1
(0, U2, Uy (5, v) (pi.j Wijo Ui.j Pi.j):
Pi..Pi.j)
(30)
where Pt.j 111. May (-11),, Uijj ~ U|-1,111 and Pt.j 1p.1111. In this dataset, P = 2 is considered, which constitutes subdomains in total. random sample of this dataset visualized in Figure D.9
25

---

t= to
t= t4
t= ta
t= t12
t= t16
1.25
1.00
0.75
3
Figure D.8: Snapshots of solution trajectories of the CE-Gauss dataset. Each row corresponds to one variable.
D.2.9 ACE
The Allen-Cahh equation [3) are considered with periodic boundary conditions and reaction coefficient 220. is initialized
1
n (x v)
Oij
sin
(31)
M
4==1
where M is random integer randomly sampled from [1,,2]], ~ 10. and asj ~ 411,11- A random sample this dataset is visualized in Figure D.10
D.2.10 Wave-Layer
The wave equation (14) is considered with the fixed-in-time propagation speed c(x, y) as piece-wise constant function in vertically layered medium with n {3, 4, 6} layers and randomly defined intersections. More details about the definition of layers are given in SM B.2.11 11 of [13]. u is initialized as a linear combination of E 2, 3, 4, 5, Gaussian distributions centered at random locations in the domain. A random sample of this dataset is visualized in Figure D.11
D.2.11 AF
In this dataset, the transonic flow is considered over different airfoils. The steady-state form of the compressible Euler equations (12) govern this system. The viscous effects are ignored and the far-field boundary conditions of Pso
1, Poc 1.0, Mr 0.8 and zero angle of attack are considered. At the airfoil boundary, no no-ppenetation conditions are considered. We refer the reader to 19 for more details about how the airfoil geometries are drawn. The flow computed on
C-grid mesh with 200 50 quadrilateral elements by second-order implicit finite volume solver. The original dataset is enhanced by adding distance field from the boundary of the airfoil which used the input of the model. The target value used in this work is the density p.
26

---

t= to
t= t4
t= tg
t= t12
t= t16
-2
0.5
Vx
0.0
-0.5
0.5
vy
0.0
-0.
0.50
0.25
Figure D.9: Snapshots of solution trajectories the CE-RP dataset. Each row corresponds one variable.
t= to
t= t4
t= t8
t= t12
t= t16
0.5
0.0
-0.5
Figure D.10: Snapshots of solution trajectories of the ACE dataset. Each row corresponds to one variable.
D.2.12 Poisson-Gauss
Poisson's s equation (6 considered unit square with homogeneous Dirichlet boundary condi- tions. the source term is composed of a superposition of M random Gaussians as
f(x. y) = Σ exp( (xr pxx,)) 2 + (v my,x)² )
M
(32)
20
where M is random integer drawn from geometric distribution, Mar.i Myyi Uj.11, and a, M10.225,0.11- A few input-output pairs of this dataset are visualized in Figure D.13
D.2.13 Elasticity
The solid body foree balance equation [77 is considered for unit square hyper-elastic incompressible Sivlin-Saunders specimen with a hole at center. The geometry of the hole randomly sampled from a certain distribution (see [19 for more details) such that the radius always takes a value between and 0.4. The specimen clamped at the bottom boundary (zero displacement) and is under constant vertical tension traction its top boundary. The stress field is computed for different
27

---

t= to
t= t4
t= t3
t= t12
t= t16
4000
-3500
Figure D.11: Snapshots of solution trajectories of the Wave-Layer dataset. Each corresponds one variable.
sample #0
sample
sample #2
sample #3
sample #4
0.8
0.6
Figure D.12: Solution samples of the AF dataset (only density). The full domain is much larger and the mesh is more sparse than the visualized frame around the airfoil.
1.5
1.0
0.5
0.0
0.005
0.000
Figure D.13: Inputs (top row, source function) and outputs (bottom row, solution) the Poisson- Gauss dataset. Each column corresponds to sample in the dataset.
geometries with finite elements solver with about 100 quadratic guadrilateral elements. The target function this dataset the stress field, which is available at 972 unstructured coordinates. The original dataset is enhanced by adding a distance field from the boundary of the hole that is given as input to the models.
sample #1 1
sample #2
sample #3
sample #4
sample #5
1e3
Figure D.14: Solution samples of the Elasticity dataset.
28

---

E
Details of the Experiments
central design principle of RIGNO is to minimize the need for extensive hyperparameter tuning by employing fixed set of hyperparameters across all tasks. To this end, we identified robust config- uration by tuning on a single internal dataset and subsequently applied the same hyperparameters to all other datasets without further adjustment. The size and performance of RIGNO are primarily influenced by only two key architectural hyperparameters: the number of processor blocks and the latent dimensionality of edge and node features. A comprehensive list of the hyperparameters used in RIGNO-18 18 and RIGNO-12 is provided in Table Since most of these hyperparameters require no task-specific tuning, the same configuration can be readily adopted for new problems or datasets. The only distinction between the two architectures lies in the number of message- message-passing blocks. Unless otherwise stated, the settings described this section are used for all experiments. As shown in Table both RIGNO- 18 and RIGNO-I 12 achieve consistently strong performance across diverse benchmarks without the need for task-specific tuning.
The last feed-forward block in the decoder that independently acts on physical nodes does not have normalization layers, i.e., it is only composed of two-layer MLP. A nonlinear activation function is applied on the output of the hidden layer in all MLPs. For the main MLPs of the feed-forward blocks, we use the swish activation function [37]; whereas for the MLPs of the layers, the sigmoid function is used.
Table E.1: Architectural hyperparameters of RIGNO and their corresponding values for RIGNO-18 and RIGNO-12.
RIGNO-18
RIGNO-12
Description
2.7M
1.9M
Total number of trainable parameters.
18
12
Number of message passing blocks.
4
4
Factor for sub-sampling the given coordinates. As an example, with 16, 000 input coordinates, get 4, 000 coordinates in the regional mesh.
160000
1.0
1.0
Overlap factor of the support sub-regions in the encoder.
2.0
2.0
Overlap factor of the support sub-regions in the decoder.
6
6
Number of regional edge levels.
128
128
Dimension of the latent node features.
128
128
Dimension of the latent edge features and mes- sages.
1
1
Number of the hidden layers in the MLP of the feed-forward blocks.
128
Dimension of the hidden layers in the MLP of the 128
feed-forward blocks.
16
16
Dimension of the hidden layers in the MLPs of conditional normalization layers.
4
4
Number of frequencies for encoding coordinates, i.e., Kiteq
0.5
0.5
Probability for edge masking.
We train our models with 1024 training samples for all datasets; except for the AF dataset where 2048 samples are used in the trainings. We use different set of 128 samples for validation, and 256 samples for testing. The test errors reported Table 2 are exceptionally calculated on 240 test samples. We train models for 2500 epochs the time-dependent datasets, and for 3000 epochs
29

---

on the time-inee -emeiiddeenddent datasets. We calculate the errors on the validation dataset regularly and keep the model parameters that give the lowest validation error. We use batch size 16, the Adam W [27] optimizer with weight decay 10-8, the mean squared error loss function, and cosine learning rate decay with linear warm-up period. The learning rate warmed up in the first 2% training epochs from 10-5 to 2 x 10-1, and decays with cosine schedule to 10-6 in the next 88% of the epochs. In the last 10% epochs, it exponentially decays to 10 For all the time-dependent datasets, the fime-derivative time stepping scheme is used due to its better performance (see Figure [c). We use the all2all training strategy [33] with a warm-up period. We first train the model with all the pairs with lead time 2ât for few epochs, then add the pairs with lead time 4åt and train for more epochs, and continue until all lead times are covered. The warm-up period for the all2all training takes 20% of the training epochs.
Most of the experiments with RIGNO have been done using NVIDIA GeForce RTX 3090 GPUs (24GB). Some experiments have been carried out with more than 2 (up to 8) GPUs in order to speed up the trainings. We estimate total of 9'000 GPU hours for reproducing the presented experiments with RIGNO (including SM).
We evaluate the estimates of the models against the ground-truth target functions following the metric used in [33.. We first normalize both the model estimate and the ground-truth function values by shifting and scaling them with global equation-wide widd mean and standard deviations. We then compute the relative L' error per each variable (group) as
Σ;
2
R.
(33)
where "i, and ug are respectively model estimate and ground-truth function values for sample m, variable group k, and spatial coordinate i. For all the Navier-Stokes and the Compressible Euler datasets, we also sum the components of the the velocity vector here. Hence, we get scalar relative error og per each target variable (group) and per each sample. We average over i to get a scalar relative error per sample, and report the median error over all samples in the validation or the test dataset.
F Further Experimental Results
There are many hyperparameters that can affect the training stability well as different aspects of the performance of trained RIGNO. the experiments presented in this work, we fix a set of architectural choices that perform fairly well for all the datasets and training sizes. These hyperparameters are listed in Table E.I Note that these hyperparameters do not correspond to the best-performing architecture, but rather they guarantee reasonable performance with all the desired properties. In Section F:3] we show that the performance can be improved by using larger models. In the experiments of this chapter,
perturb or few of these hyperparameters while leaving the rest fixed in order to study how the performance of the model changes with respect to the changed hyperparameters. The trainings carried out using the settings described in Section E unless otherwise stated.
F11 Training randomness
In order to quantify the dependence of the final results on training randomness, we train RIGNO five times, each time with a different seed for the pseudo-random number generator. We store the model with the lowest validation error during training as in all experiments, and check for fluctuations in the test error across all five repetitions of each experiment. In these experiments, we use the unstructured version of the datasets (random point clouds), we train using 512 training samples with batch size 8, and do not fine tune with fractional pairing. The results are summarized in Table F. where the maximum and minimum test errors across the repetitions are reported for each dataset. These fluctuations are attributed to the randomness processes such as weight initializations, dataset shuffling, edge masking, and random sub- aub-sampling of the coordinates. In the rest of the experiments, we only train each model once with a fixed training random seed, keeping in mind that the errors can slightly deviate because of training randomness.
30

---

Table F.I: Statistics of relative L test errors with different random seeds. We train RIGNO-12 on 5 512 samples from each datasets. Each experiment repeated 5 times.
Dataset
Error [%]
Mean + Standard deviation
NS-PwC NS-SVS
2.50 + 0.157 1.05 = 0.040
F.2 Number of regional nodes
The number of regional nodes is independent from the given physical nodes and can be chosen arbitrarily. With more regional nodes, the support sub-regions will be smaller on average, leading to fewer edges in the encoder and the decoder. However, the processor steps become significantly more expensive. In order to investigate the effect of different numbers of regional nodes on the performance, we train slightly modified RIGNO-18 18 on the NS-SI. and the NS-SVS datasets. The test results are presented Table F2 With only 1024 regional nodes, the model already achieves reasonable error for both datasets. With more regional nodes, the error follows decreasing trend with both datasets, although less significant with NS-SVS. However, the inference times show that fewer regional nodes are favorable in terms of computational efficiency. RIGNO-18 18 is 5 times more time-efficient with 1024 regional nodes as compared to 16384 regional nodes. Moreover, the decreasing trend of the error breaks with too many regional nodes. With sub-sample factor 1.0 (16384 regional nodes), the support sub-regions become so small that most of them only contain one physical node and hence cannot reflect good representation of the whole spatial region. In the main experiments of the current work, we have chosen a sub-sample factor of 4.0 to maintain a balance between accuracy and efficiency.
Table F.2: Test errors and inference times with different number of regional nodes. The experiments are conducted with a variation of RIGNO-18, trained 256 training trajectories for 2000 epochs without all2all pairing. Although the sub-sampling factor for obtaining the regional nodes changes, multi-scale edges are always obtained based on sub-meshes of factor 2.0 until the regional mesh saturates.
Regional nodes
Median relative L error [%]
Inference time [ms]
(sub-sample factor)
NS-SL
NS-SVS
1024 (16.0)
4.32
2.90
7.0
2048 (8.0)
3.96
2.93
9.0
4096 (4.0)
3.49
2.84
12.9
8192 (2.0)
3.33
2.83
28.3
16384 (1.0)
4.12
3.05
35.8
F.3 Scaling with model size
In order to evaluate how the performance of RIGNO scales with model size, we start with RIGNO-12 as the base architecture and change the latent sizes (all at the same time), the number of message passing steps in the processor, the overlap factors (both at the same time), and the number of times that the regional mesh is sub-sampled when processor edges are built (number of regional edge levels). 256 samples of the PS-PwC dataset without fractional pairing and calculate
We train these models the error at 11 on 128 test samples by using the AR-2 rollout scheme (defined in Section results are plotted in Figure
The
We observe that the latent sizes and number of processor steps have the largest effect on the model capacity. The error consistently falls down when these hyperparameters are scaled, and we expect to get even lower emoor with larger RIGNO architectures. With larger overlap factors, although the number of edges (and hence the computational cost) increases rapidly, the rate of convergence is relatively low. Therefore, the authors recommend using the minimum overlap, especially for the encoder. With more regional mesh levels, the other hand, the number of edges only increases
31

---

logarithmically, which is also reflected in the inference times. With 4 regional mesh levels, we obtain 3.75% error, which is extremely lower than the 17.10% error with only a single regional mesh level. The inference time, however, is only increased by 14%. These results show the effectiveness of longer edges the construction of the regional mesh. However, depending on the space resolution of the inputs, there is saturation limit for this hyperparameter. If we consider a downsampling factor of 4, with 16'000 input coordinates, for instance, we end up with a mesh with only 4 nodes after downsampling steps. With so few nodes we can only get 6 additional edges at most, which explains why the decreasing trend the error stops at 4 downsampling steps.
Inference time
fevvnne time
Inference time (us)
Inference [uss] 10.0
4.12
10.2
7.57
18.
10.2
20
Benor (36)
Bevor
Frvor
3.0
2.0
Latent ston
Proceparr stepe
Overlop foctor
Regional mosh lnvele
Figure F.1: Test errors at / = 41 with different strategies for scaling model size. The experiments are carried out on the NS-PwC dataset.
F.4 Spatial continuity: resolution invariance
In Figure 3a we presented results of zero-shot super-resolution and zero-shot sub-resolution infer- ence of RIGNO. The results show that RIGNO maintains its accuracy when inferred with higher resolutions but the accuracy drops with lower resolutions. When trained with edge masking, zero-shot generalization to lower resolutions is considerably improved. Generalization to higher resolutions often desirable feature since the model can be trained with lower resolutions with considerably less training costs, and be deployed for high-resolution inputs and outputs. In Figure F22 we present the results of the resolution invariance test MeshGraphNet. The results show that MeshGraphNet admits great dependence on the training resolution, with very large errors on any resolution except for the training resolution. This dependence is alleviated when MeshGraphNet is trained with edge masking. The accuracy is much less affected both with lower and higher resolutions. However, edge masking alone is not enough make MeshGraphNet resolution invariant. With only 4 times higher resolutions, the error is almost doubled 16%) as compared to the error with the training resolution (30.36%)
100
Error [%]]
30
40
FM
50% EM
42
522
64
76
922
128
Inference resolution
Figure F.2: Resolution invariance test if MeshGraphNet with 642 training resolution on the NS-PwC dataset. EM stands for edge masking. The y-axis shows the relative test error (at t = (14) when the spatial resolution of the inputs and outputs varies from the training resolution.
F.5 Temporal continuity: generalization to unseen time steps
In order to assess the generalization of the trained models to unseen lead times т. we train RIGNO on the NS-PwC dataset with 2A/ time resolution, and test it the test dataset with At time resolution. With the incorporated mechanisms, namely time-conditioned normalization layers, all2all training, and fractional pairing fine-tuning, expect to see monotone increase in the error as we increase
32

---

the lead time; which is very well captured in our tests as visualized in Figure F3 Here, the lead time is fixed at each line in the plot, and the imput time snapshot is changed. This explains why the first entries are missing for larger lead times. In the uppermost line (r = 921), for instance, the first entry corresponds to 0, = 9.1; the second entry to 1 т = 10 and so on.
In Figure F.3 the errors increase with larger lead times, which is expected since estimating the solution of the PDE with a higher lead time is inherently more difficult. In these experiments, we have removed the pairs with 7 > 8At from the all2all training pairs. Hence, T 9 is considered as extrapolation in lead time. Since lead times r = 3åt, т 5 and T 7åt are also not present in the all2all pairs, they are considered as interpolation in lead time. We can see that the monotonicity of the increase in the error is preserved for the extrapolated as well as the interpolated lead times. Moreover, the error increases in controlled manner. The authors believe that this property attributed to incorporation of condiiontt normalization layers in the architecture and the use of time-derivative time-stepping strategy. However, with the smallest lead time r = At, it was observed that the monotonicity was not preserved without fine-tuning with the fractional pairing strategy (see Figure 3b).
3.5
IA
r
2At
3.0
3A
4A
2.5
SAt
Enror [%]]
6A
2.0
7A
8A
1.5
T 9A
1.0
trained
interpolaced
0.5
extrapolated
0.0
0
2
6
10
12
14
16
18
20
Time index
Figure F.3: Direct test errors with different 7 for the NS-PwC dataset. The vertical line indicates t = f14, which is the last time snapshot seen by the model during training. Different lines denote different lead times. The gray lines correspond to lead times present in the training dataset; the green lines highlight interpolation in lead time; and the red lines indicate extrapolation in lead time. The model is trained 512 samples with batch size 4.
F.6
Model uncertainty
The edge masking technique randomly disables subset of the edges in each cesssag-passsing block. By adding randomness to the model, edge masking allows inferring the same model multiple times with the same input in order to generate an ensemble of estimates at a later time. The standard deviation of the generated ensemble can be interpreted as an indicator of model uncertainty on each output coordinate.
Figure F.4 visualizes the mean and the standard deviation of an ensemble of estimates at f = fu generated from input at t to. along with the absolute error of the mean with respect to the ground- truth solution. More such visualizations are available in Section H These visualizations show how the regions with the highest errors largely overlap with the regions where the ensemble has the highest deviations; which supports the interpretation of the deviations as a measure for model uncertainty, and shows that they good indicators for inaccuracies in the estimates generated by the model.
Generally, the uncertainty of the model increases over autoregressive rollouts or with larger lead times. The model uncertainty can therefore be used to optimize the set of lead times in an autoregressive path. A possible algorithm would be to adaptively choose the largest lead time 7 that results in uncertainties smaller than threshold. Since model inference is fast and differentiable, grid-search
33

---

algorithms can also be implemented. Similarly, the input discretization can adaptively be refined at each time step in the areas where the uncertainty of the model in its output is the highest.
Ensemble mean
Absolute eror
Ensemble std
Figure F4: Mean and standard deviation of an ensemble of autoregressive estimates at 114 for test sample of the NS-SL dataset.
E.7 Robustness to noise
In many applications, a learned operator needs to inferred with inputs that are collected experimen- tally. Since measurement noise is inevitable in any experimental setup, robustness noisy inputs is vital for such solution operators. Ideally, the operator should dampen the noise in the input; or at least the noise should not propagate in uncontrolled way. To evaluate the robustness of RIGNO input noise, we train it using clean, noise-free data, infer it with noisy inputs, and compute the error of the output with respect to clean ground-truth data. This setting emulates applications where the operator has been trained on noise-free free high-fidelity simulation data and is applied on noisy experimental data.
noise O impus
Noise-ind ind-ced error [%)]
109
Noiso-inducell error [9]]
10a
0.5
2.0
10-1
10-1
Tuoke Imput
0.5
2.0
10-2
10
0
2
6
10
12
14
2
4
8
10
12
14
Time index
Time index
(a) Heat-L-Sines
(b) NS-PwC
Figure F.5: Autoregressive errors introduced by the noise the input (relative). The y-axes show the difference in the relative error with noisy and clean inputs. The solid lines (best) show the results with edge masking and overlap factors 0. The dashed lines show the results with the same overlap factors but without edge masking during training. The dotted lines show the results with edge masking and overlap factors 1. All models are trained on 256 trajectories, without all2all training. and without fractional pairing fine-tuning
34

---

For each input in the test dataset, we compute its standard deviation O'pput. and add a Gaussian noise with standard deviation noise inpust with noise level Cnoise {0.005, 0.02). Starting with a noisy input, we apply the operators autoregressively with Tmax 2 The errors in the resulting trajectories are partially due to the input noise, and partially to the approximation error of the operator. In order to distinguish between these two, we also unroll the operator with clean input (Cnoise = 0), and subtract the resulting errors from the errors with noisy inputs. This procedure gives a measure for the part of the error that is introduced by the noise in the input, denoted by noise-induced error. We present the relative noise induced error with respect to the approximation error: |Cnoisy Cclesal/ / /elgann where represents the autoregressive test error at t f14- The results are visualized in Figure|F.3 With 2% input noise at tu, the noise is very well dampened to only 0.5% at t14 of the Heat-L-Sines dataset. With larger support sub-regions, the regional nodes aggregate data from a larger pool of physical nodes. Therefore, they are less sensitive to noisy inputs, and a stronger noise dampening is observed. Edge masking also plays important role in noise dampening. When larger support sub-regions are used along with edge masking, 2% input noise is dampened below 0.2% at 114-
Furthermore, noise dampening of RIGNO is not limited to Gaussian noise, as there has been no indication of this type of noise during training of the model. The authors believe that RIGNO is also robust to the errors that are introduced during unrolling, if we consider the introduced errors as noise with certain distribution. In other words, RIGNO corrects itself within successive autoregressive steps. Therefore, for certain problems, small -but not too small- lead times might be advantageous to larger ones due to noise dampening.
F88 Out-of-distribution generalization
central challenge in operator learning, and in supervised machine learning more broadly, is achieving strong generalization to out-of-distribution (OOD) inputs and outputs. These models often experience a degradation in performance when evaluated on samples that deviate from the training distribution. While some degree of performance decay is inevitable, robust framework should maintain reasonable accuracy under mild to moderate distribution shifts. To evaluate RIGNO's robustness in this regard, we train it on three datasets taken from 38). namely Smooth Transport, Darcy Flow, and Compressible Euler Equations. The trained models are then tested on 128 in- distribution (ID) and 128 OOD samples. The OOD samples differ from the training distribution in the location of perturbations, characteristic length scales, or geometry complexity, All models are trained for at least 2000 epochs using a batch size of 4 to ensure convergence, and all evaluation metrics are computed on unnormalized (raw) test data.
The results are summarized in Table F3 Interestingly, RIGNO exhibits strong robustness to OOD samples in the Darcy Flow and Compressible Euler problems, but experiences pronounced perfor- mance drop for the Smooth Transport problem when evaluated on OOD inputs. In the case of the Compressible Fuler problem, the OOD samples feature airfoils with more complex geometries, while in Darcy Flow, the OOD inputs differ primarily in their characteristic length scales. Both cases can be interpreted as involving local variations in the structure of the input and output functions, to which RIGNO generalizes well. In contrast, the Smooth Transport problem introduces non-local shift, as the center of the Gaussian perturbation changes significantly in the OOD samples, leading to the observed performance degradation.
Table F.3: Median relative L errors [%] on 128 in-distribution (ID) and 128 out-of-distribution (OOD) test samples. The errors are calculated on unnormalized samples.
Dataset
Training samples
ID
OOD
Darcy Flow
256
2.57
3.18
Smooth Transport
512
0.72
6.50
Compressible Fuler
512
0.51
0.62
F.9 Pre-training and transfer learning
PDE foundation models have been proposed 13
30 as a solution to overcome the high data
demands of neural operators. In this context, we have assessed the generalization capabilities of
35

---

RIGNO to new input and output distributions. To this end, we train RIGNO on 1024 trajectories (2000 epochs, without all2all pairing) of one dataset (e.g., NS- Gauss) and transfer dataset (e.g., NS-PwC). We fine-tune the pre-trained RIGNO with only few trajectories of the target dataset and compare the results with a randomly initialized RIGNO. We have chosen the NS-Gauss, NS-PwC, and NS-SL datasets, which have the same underlying PDE (incompressible Navier-Stokes equation) but correspond to very different data distributions (see Figures D.4 D.5 and [.6]. The results are depicted in Figure [6 For the NS-PwC dataset, we find that with only 4 trajectories, the transferred RIGNO pre-trained on NS-Gauss has an error of 13. %, while a RIGNO trained from scratch requires 128 trajectories achieve the same error, providing 32x gain in sample efficiency with transfer learning. Notably, the gain is smaller (roughly 4x) when the model is pre-trained on the NS- SL dataset. Similarly, for the NS-SL dataset as the target dataset, transfer learning with RIGNO
different
32 samples lead to an error of 4.7%, which is 7 times less than the trained-from-seranne (33.4%) with 32 samples. These examples illustrate the ability of RIGNO to transfer very well out of distribution.
Initialization
Initialization
random
NS-Gauss
NS-SL
random
NS-Gauss
NS-PwC
100
100
50
50
Error [96]
Error [9]]
2
2
1
128 Number of training trajectories
2
61
256
16 32 Number training trajoctories
256
(a) Target dataset: NS-PwC.
(b) Target dataset: NS-SL.
Figure F.6: The effect of pre-training on data efficiency of RIGNO-12 with the unstructured NS datasets. All errors are relative test errors on 256 test samples.
F.10 Performance benchmarks
We benchmarked RIGNO against the best-performing baseline, GINO, in terms of training time, accuracy, inference time, and peak memory usage during inference. For training, we used the maximum feasible batch size for RIGNO and fixed batch size of 64 for GINO. In cases involving varying geometries across training samples, such as in the Elasticity dataset, the original GINO implementation does not support batch processing, requiring reduced batch size of 1. We evaluated two RIGNO-18 configurations with 96 and 128 latent features, respectively, and two GINO variants comprising and FNO layers.
Table F.4: Performance benchmarks on two unstructured datasets evaluated on a single NVIDIA GeForce RTX 3090 GPU (24GB). TT: training time [s]: TE: test error [%]; IT: inference time [ms]; and PM: peak memory usage [MB].
Miisoon-Gauss
Elasticity
Architecture
Size
TT
TE
IT
PM
TT
TE
IT
PM
GINO (3)
3.6M
1.33
5.24
7.8
370
18.2
3.36
4.1
80
GINO (7)
8.4M
1.72
8.01
10.1
406
27.1
3.94
6.4
116
RIGNO-18 (96)
1.4M
24.9
2.23
8.1
121
2.4
4.33
1.5
51
RIGNO-18 (128)
2.5M
28.4
2.09
12.7
159
3.1
3.71
1.9
64
The results, summarized in Table F.4.
indicate that RIGNO trains more slowly than GINO on datasets with fixed geometries and input coordinates, but trains faster on datasets with variable geometries, primarily due its support for batch processing. Although longer training times are drawback.
36

---

training constitutes one-time cost, whereas inference speed and memory usage are typically more critical in practice. In this regard, RIGNO and GINO exhibit comparable inference times, yet RIGNO consistently achieves substantially higher accuracy. Furthermore, RIGNO is up to three times more memory-efficient than GINO. Owing its multi-scale architecture and computational efficiency, RIGNO stands out as powerful end-to-end GNN that delivers state-of of-e-of-the-art accuracy and robustness, while remaining competitive with strong non-GNN baselines in both speed and memory consumption.
G Baseline Architectures
G.1
CNO
A Convolutional Neural Operator (CNO) is model that (approximately) maps bandlimited functions to bandlimited functions (38. Let B., be the space of bandlimited functions with the bandlimit w. A
CNO is compositional mapping between function spaces G B. (D) + (w( and is defined as
g: P(u)
00
+ Q(v..) u.
(34)
where
U++ P: E K((o))
-
(35)
where I. is the number of CNO blocks and D is the domain.
Here, P is the lifting operator, Q is the projection operator, either the upsampling down- sampling operator, Kr is the convolution operator and Σ₁ is the (modified) activation operator. For continuous time conditioning, we used lead tnne-condtiooeed instance normalization, defined by
(ee (x) a(t) (N)v B(t)
(36)
where input function, (v) is regular instance normalization, and arl)) and 3(1) are (small) MLPs.
In all the benchmarks, the CNO architectures and training details utilized in this work are taken from the paper [13] and are as follows:
Trainable parameters: 39.1M
Number of epochs: 400
Lifting dimension:
54
Up/downsampling layers:
Residual blocks the
bottleneck:
Residual blocks in the middle layers:
Batch size:
32
Optimizer:
AdamW
Weight decay:
10-6
Initial learning rate:
5 10-4
Scheduler: Linear with decreasing factor of 0.9 every 10 epochs
Early stopping:
If the validation loss does not improve for 40 epochs
G.2
ScOT
The scalable Operator Transformer (scOT) is a multiscale transformer model mapping discretized functions regular cartesian grid to discretized functions [13]. It patches the input functions and
37

---

projects these visual tokens linearly into latent space. A Swin [25 encoder-decoder structure with ConvNeXt [26] residual blocks transforms the latent tokens in hierarchical manner. SwinV2 transformer blocks apply cosine attention not on the full token grid, but only on windows of tokens, effectively reducing the computational burden of full attention. Downsampling operations in the encoder after each level merge four tokens into one, therefore coarsen the computational grid; upsampling in the decoder reverses these operations. After the decoder, the latent tokens are mapped back linearly into the physical space.
The model accepts continuous time or lead time inputs. is conditioned through every single layer norm module by that lead time. The conditioned layer norm N applied to latent token v E IRC at some level is implemented as follows:
m(v) 1
Eiv (v C
((v) ((v)--((v)
(37)
1=1
u(v) 10
-(v,t)=al) = a(l) 0 d²(v)
3(t)
(38)
with lead time t Rzo, 1e E R vector of ones, and (t)), ((t) being learnable affine maps from lead time to gain and bias.
As with CNO, we utilize the model and training setup from 13 For the sake of completeness, they are given as follows:
Trainable parameters: 40M
Patch size:
Embedding/latent dimension: 48
Window size: 16
Number of levels:
Number of SwinV2 blocks each level:
Number of attention heads in each level:
63,6,
Number of ConvNeXt blocks in the residual connection:
2
Optimizer:
AdamW
Scheduler:
Cosine Decay with linear warmup of 20 epochs
Maximum learning rate:
10-4
Weight decay:
10
Batch size:
40
Number of epochs: 400
Early stopping: If the validation loss does not improve for 40 epochs
Gradient clipping (maximal norm):
G.3 FNO
The Fourier neural operator (FNO) is a model that efficiently maps bandlimited functions to bandlim- ited functions convolutions Fourier space 20 and the fast Fourier transform (FFT). Similar to CNO, the FNO is composition of lifting, convolution, and projection operators. The convolution operator, KG (o), may be viewed as multiplication of data, U, learnable kernel function, Rø. in Fourier space. Here, o learnable parameters of neural network. The convolution operator
38

---

may therefore be expressed as,
(K(()))) F= F- (Rs (Fo,)) (z),
(39)
where
11+1 (z) (a IN) (We((a) (K(()))) (z)), -SI5
(40)
W is a local integral operator and F is the Fourier transform, IN is (lead condicondii instance normalization. For efficiency, the Fourier transform is implemented via the fast Fourier transform (FFT) with a computational complexity of Oog n). where n is the number of points in the discretized mesh. Likewise, low-pass filter removes high-frequency components of the Fourier transform, such that the weight tensor Rg contains only kmar < n. modes. Because kmaz is independent of n, generally taken as a hyperparameter between and 32, the inner multiplication has a complexity of O(kmaz). As with CNO, we condition the model on lead time through instance normalizations.
The same setup as in 16 is used for this work:
Trainable parameters: 37M
Lifting dimension: 96
Number of Fourier layers:
Number of Fourier modes 20
(kmax):
Optimizer:
AdamW
Scheduler:
Cosine Decay
Initial learning rate:
10-4
Weight decay:
10-6
Number of epochs: 400
Batch size:
40
Early stopping: If the validation loss does not improve for 40 epochs
G.4 Geo-FNO and FNO DSE
The FFT is limited to equispaced rectangular grids, and this restriction is inherited by the FNO. To apply the FNO to general geometries, the geomelry-aware Fourier neural operator (Geo-FNO) [99 learns diffcomorphic deformations between general point clouds and such rectangular domains. Let Dg be the physical domain, Dc the computational domain, I E Da. and E DE be the corresponding mesh points, Ta = (x) C Da the input meshes, T = {g(()) C D the computational meshes. The diffeomorphism d. transforms the points from the computational mesh the physical mesh by
Qa : D + Da.
(41)
{ + I,
with inverse o- and both oa and --1 are smooth. The Fourier transform is subsequently approxi- mated as
(Fau(k))
Σ ( (x).A) TET
(42)
T
where m(x) is weight matrix based the distribution of points from the input mesh.
In contrast, the FNO with discrete spectral evaluations (FNO DSE) 24 formulates a band-limited discrete Fourier transform based on the input mesh. Based the observation that the number of Fourier modes employed by the FNO, kmazs is small and independent of the number of mesh points, direct matrix-vect multiplication can be performed with (kkooz)) complexity. For small
39

---

Kmaz, this operation is more efficient than FFT. The discrete Fourier transform matrix for mesh points Ta = {x(h)} C. Da is constructed by
VAJ 1 Vn
kmax- 1,n-1
(43)
R.==0
Using this matrix as a basis for the bandlimited Fourier transform, the FNO DSE avoids the construc- tion of computational mesh. Instead, it works directly the physical mesh.
In all the benchmarks, the Geo-FNO and FNO DSE architectures and model hyperparameters are taken from 24 and fixed for all experiments follows:
Trainable parameters:
3.68M (Geo-FNO), 3.30 (FNO DSE)
Number of epochs:
1000
Lifting dimension:
48
Number of Fourier modes
14
(kmar))
Number of Fourier layers:
Computational mesh size
40 x 40
(Geo-FNO):
Batch size:
16
Optimizer:
Adam
Weight decay:
10-5
Initial learning rate:
10-4
Scheduler:
Linear with decreasing factor of 0.9 every 10 epochs
Early stopping:
If the validation loss does not improve for 40 epochs
Based on the works [24 19]. two sets of favorable training hyperparameters were tested for each
experiment. The results of the best model are reported. The first set of hyperparameters is:
Initial learning rate:
1 10-3
Scheduler: Linear with decreasing factor of 0.5 every 50 epochs
The second set is:
Initial learning rate:
5 10-3
Scheduler: Linear with decreasing factor of 0.97 every 10 epochs
G.5 GINO
Geometry-Informed Neural Operator (GINO)
23 composed of three main operators:
G
E P D.
(44)
where E is the encoder, P is the processor, and D is the decoder. Specifically, & maps the input data from an unstructured mesh to a uniform latent grid using graph neural operator (GNO). P processes the latent-gria features via the FNO. D maps the latent features back to the desired output mesh (or query points) using GNO.
GNO layers approximate an integral operator on a local neighborhood graph. Given a ball of radius r > 0 around each point D,
07(x) = ^(x.8) 11)( dy.
(45)
B,))))
40

---

we discretize the integral by selecting M neighbors, {yi.... , C (x)x), and form a Riemann sum:
M
(x)
Σ Adr,yi) 1y-) ((3s).
(46)
where vu-) represents the local weights for integration. These localized integral updates are performed across a graph that connects each node z to its neighbors in (z)..
In all the benchmarks, the GINO architectures and model hyperparameters are taken from [X] and
fixed for all experiments as follows:
Trainable parameters: 3.7M
Number of epochs:
500
Radius for GNO layers:
0.033
GNO transform type:
Linear
Projection channels:
256
Lifting channels: 256
FNO input channels: 64
Number Fourier of modes: 16
Number of Fourier layers: 3
Optimizer:
Adam
Batch size:
64
Initial learning rate:
1- 10
Scheduler: Linear decay with factor of 0.8 every 100 epochs
G.6 MeshGraphNet
The MeshGraphNet 33 architecture is an end-to-end graph-based neural operator that performs message passing directly on the given coordinates (physical nodes). Since the architecture shares many similarities to RIGNO, we can obtain MeshGraphNet by making few modifications to our architecture. We consider RIGNO-12 and perform the downscaling (encoder) and upscaling (decoder) processes by considering very small overlap factors, i.e., very small support radii. Regional nodes lose their meaning since each of them is only connected to one corresponding physical node. Removing the multi-scale edges in this architecture leads to a MeshGraphNet with 1.9M parameters. We train this MeshGraphNet with edge masking to be more consistent with RIGNO but use larger batch size (64) to compensate for the other hyperparameters and settings are the same as reported in Section
training time due removing the regional nodes. All
G.7 UPT
The Universal Physics Transformer (UPT) is a unified neural operator framework designed for flexible and scalable learning across wide range of spatio-temporal problems, including both Eulerian (grid- based) and Lagrangian (particle-based) simulations [I_ UPT consists of three components: an encoder that compresses the input into a latent representation, a transformer basedd approximator that propagates dynamics entirely in the latent space, and decoder that queries this representation at arbitrary spatial locations. The model is conditioned on time and boundary conditions via DiT-style feature modulation.
We adopt lightweight configuration of UPT in our experiments. Specifically, we use 64 latent tokens and 64 hidden dimensions. For time-dependent datasets, this results in 2. 18M trainable parameters, and for time-independent datasets (without time conditioning). the number parameters 0.74M. A larger variant with 256 latent tokens and 192 hidden dimensions with 10.7M parameters has also been considered. We found that the large UPT strongly suffers from optimization difficulties on the datascts considered in this work, leading to higher errors than the ones reported in Table with all
41

---

datasets. In all the presented benchmarks, UPT architecture and model hyperparameters are fixed for all experiments as follows:
Trainable parameters:
2.18M or 0.74M (see above)
Number of epochs:
500
Number of latent tokens:
64
Hidden dimensions:
64
Number of supernodes: 2048
Encoder and decoder depth:
4 layers each
Approximator depth:
4 layers
Number of attention heads:
Optimizer:
AdamW
I Learning rate:
1.10-3
Scheduler:
Cosine decay
Batch size:
32
Early stopping:
If validation loss does not improve for 40 epochs
G.8 Transolver
Transolver [45 is an attention- neural operator framework that learns global latent tokens from input features and processes them through self-attention mechanisms. In our experiments, we observed that this architecture scales poorly to large-scale datasets. This limitation arises from its of global interactions implemented as independent MLPs between input points and sparse tokens. Consequently, very large number of MLP evaluations are required, quickly exhausting GPU memory. For instance, on an RTX 4090 GPU, the maximum feasible Transolver model size during training is approximately 0.5M parameters. For time- lime-independent datasets, we were able to increase the model size to 0.9M parameters, The corresponding results are summarized in Table G. While both RIGNO variants exhibit clear advantages on time-dependent problems, the accuracy is comparable for elliptic problems.
Table G.I: RIGNO and Transolver benchmarks on datasets with unstructured space discretizations. All training datasets have 1024 trajectories except for the AF dataset (2048).
Dataset
Median relative 7,1 error [%]
RIGNO-18
RIGNO-12
Transolver
Wave-C-Sines
5,35
6,25
27.7
CE-RP
3.98
4.92
17.8
Poisson-Gaute Gauss
2.26
2.52
2.02
AF
1.00
1.09
1.21
Elasticity
4.31
4.63
4.92
H
Visualizations
Visual inspection allows getting more insights about several aspects of the performance of RIGNO. Estimates and statistics of ensembles of estimates are visualized in this section for all datasets. All architecture hyperparameters and training details correspond to Section All the presented samples are from the test datasets, and have been obtained by applying the operator autoregressively with Tmax 2A/ (where applicable).
Model estimates and their absolute errors, along with model input and the ground-truth output are visualized in Figures
H.I H.21 Visual inspection of model estimates for the CE-Gauss and CE-RP
42

---

datasets reveals that there are two main reasons that explain the high errors for these datasets. First, model estimates are generally smoother than the ground-truth solutions. This inability produce high-frequency estimates is reflected in the high errors where vortices are present. This is evident in the high-frequency details of velocity components in these model estimates. Second, although the model is able to capture shocks quite well at the right locations, it generates smoother transition with smaller gradients. This is very well witnessed the absolute estimate errors of Figure H.19 where the highest errors are located in small neighborhoods around shocks. The statistics of ensembles of estimates are visualized Figures H.22 H.25 Estimate ensembles are generated by unrolling the
operator randomly 20 times. Because of random edge masking, the estimates will slightly be different every time. The statistics of the estimate ensembles can act as indicators for model uncertainty. See Section for a broader discussion.
Input
Ground-truth
Model estimate
Absolute error
6
0.0
2.5
5.0
7.5
Figure H.l: Estimation RIGNO-18 for a random test sample of the Heat-L-Sines dataset at f14.
Input
Ground-truth
Model estimate
Absolute error
-5
2
0.5
1.0
1e-2
1e-2
1e-2
Figure H.2: Estimation of RIGNO-18 for random test sample of the Wave-C-Sines dataset at t14-
43

---

Input
Ground-truth
Model estimate
Absolute error
-1.0
-0.5
0.0
0.5
0.0
1.0
1e-1
-1.0
le-1
Figure H.3: Estimation of RIGNO-18 for random test sample of the unstructured NS-Gause dataset
at /14.
Input
Ground-truth
Model estimate
Absolute error
6
-2
0.5
1.0
16-1
18-1
1e-2
-5
0.0
0.5
1.0
1.5
16-1
1e-1
1e-2
Figure H.4: Estimation of RIGNO-18 18 for a random test sample of the unstructured NS-PwC dataset /14.
44

---

Input
Ground-truth
Model estimate
Absolute error
1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
1
-0.5
1.0
1e-1
1e-2
Figure H.5: Estimation of RIGNO-18 for random test sample of the unstructured NS-SL dataset
414.
Input
Ground-truth
Model estimate
Absolute error
-2
2
3
2
2
2
le-1
1e-1
1e-3
1e-1
16-1
1e-3
Figure H.6: Estimation of RIGNO-18 for random sample of the unstructured NS-SVS dataset
at (14-
45

---

Input
Ground-truth
Model estimate
Absolute error
0.9
1.0
1.1
0.8
0.9
1.1
1.
14
1.5
0.0
0.5
1.0
1.5
1e-1
-15
-1.0
-0.5
0.0
0.5
1.0
1.5
1e-1
1e-1
4
2.0
2.5
3.0
3.5
4.0
4.5
6
lc-7+2.5
1e-1
Figure H.7: Estimation of RIGNO-18 18 for random test sample of the unstructured CE-Gaus dataset 114-
46

---

Input
Ground-truth
Model estimate
Absolute error
1.3
14
1.0
1.5
2.0
3.0
3.5
4.0
1.0
-2.5
0.0
2.5
1.000
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
1e-1
1e-1
1e-1
1e-1
1.5
2.0
2
2
1e-1
Figure H.8: Estimation of RIGNO-15 18 for random test sample of the unstructured CE-RP dataset at t14-
47

---

Input
Ground-truth
Model estimate
Absolute error
775
-5.0
-2.5
0.0
2.5
5.0
7.5
1e-3
1e-4
Figure H.9: Estimation of RIGNO-12 for a random test sample of the unstructured Poisson-Gauss dataset.
Input
Ground-truth
Model estimate
Absolute error
-5
77.5
-5.0
-2.5
0.0
2.5
5.0
7.5
1e-1
1e-1
1e-3
Figure H.10: Estimation of RIGNO-18 18 for random test sample of the unstructured ACE dataset t14-
Input
Ground-truth
Model estimate
Absolute error
le-1
le-2
Figure H.11: Estimation of RIGNO-18 for random test sample of the unstructured Wave-Layer dataset at f14.
Input
Ground-truth
Model estimate
Absolute error
0.0
0.5
1.0
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
0.0
1.0
Figure H.12: Estimation of RIGNO-18 for random test sample of the AF dataset.
48

---

Input
Ground-truth
Model estimate
Absolute error
2
4
8
2
1e-1
1e2
le1
Figure H.13; Estimation of RIGNN--88 for random test sample of the Elasticity dataset.
Input
Ground-truth
Model estimate
Absolute error
-1.0
-0.5
0.0
0.5
1.0
2
1e-2
1e-2
Figure H.14: Estimation of RIGNO-18 for random sample of the NS-Gauss dataset (14-
49

---

Input
Ground-truth
Model estimate
Absolute error
-6
-
0
2
0.0
0.5
1.0
1,5
-5
0
0.0
2.5
5.0
7.5
1e-1
le-1
1e-3
Figure H.15: Estimation of RIGNO-18 for a random test sample of the NS-PwC dataset at f14.
Input
Ground-truth
Model estimate
Absolute error
1155
-1.0
-0.5
0.0
0.5
1.0
1.5
2
-5
0
-1.0
-0.5
0.0
0.5
1.0
2
1e-1
1e-2
Figure H.16: Estimation of RIGNO-18 18 for a random test sample of the NS-SL dataset at t
50

---

Input
Ground-truth
Model estimate
Absolute error
-2
2
4
3
2
2
2
4
1e-3
16-1
16-1
16-3
Figure H.17: Estimation of RIGNO-18 18 for a random test sample of the NS-SVS dataset at 414-
51

---

Input
Ground-truth
Model estimate
Absolute error
1e-2
1155
-1.0
-0.5
0.0
0.5
1.0
1.5
2
1e-1
-7.5
-5.0
-2.5
0.0
2.5
7.5
0.5
1.0
1.5
1e-1
1e-1
2.50
2.75
3.25
3.50
3.75
le-742.5
Ie-1
Figure H.18: Estimation of RIGNO-18 for a random test sample of the CE-Gauss dataset at f14-
52

---

Input
Ground-truth
Model estimate
Absolute error
1.3
1.4
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
0.0
2.5
5.0
7.5
1e-1
0.0
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
1e-1
1e-1
-2.5
0.0
2.5
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
1.5
2.0
0.0
0.5
1.0
1e-1
1e-1
1e-1
Figure H.19: Estimation of RIGNO-18 for a random test sample of the CE-RP dataset (14-
Input
Ground-truth
Model estimate
Absolute error
-5
77.5
-5.0
-2.5
0.0
2.5
5.0
7.5
0.5
1.0
1e-1
1e-1
le-3
Figure H.20: Estimation of RIGNO-18 for random test sample of the ACE dataset at f14-
53

---

Input
Ground-truth
Model estimate
Absolute error
2
16-1
1e-2
Figure H.21: Estimation of RIGNO-18 for random test sample the Wave-Layer dataset 144-
Ensemble
Absolute
Ensemble std
-0.5
0.75
1.00
1.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
10-1
Figure H.22: Model uncertainty with RIGNO-18 for random test sample of the NS-Gauss dataset.
Ensemble mean
Absolute error
Ensemble std
0.6
0.8
1.0
1.2
0.0
0.5
1.0
1.5
1e-2
16-2
Figure H.23: Model uncertainty with RIGNO-18 for random test sample the AF dataset.
54

---

Ensemble mean
Absolute error
Ensemble
0.6
0.8
1.0
1.2
0.0
1.0
15
1e-1
1--1
0.0
0.5
1.0
1.5
2.0
10-1
16-1
1e-1
Figure H.24: Model uncertainty with RIGNO-18 for a random test sample of the CE-Gauss dataset.
55

---

Ensemble mean
Absolute error
Ensemble std
1.5
0.50
0.75
0.00
0.25
0.50
0.75
1.00
Figure H.25: Model uncertainty with RIGNO-18 for random test sample of the CE-RP dataset.
56

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The relevant experimental results can be found in Tables 12 (accuracy on
multi-scale problems), Figure 3 (spatio-temporal resolution invariance), SM Section R2
(scalability), and SM Section F.4 (resolution invariance).
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The limitations are discussed in the last paragraph of the Discussion section.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
57

---

Answer:
[NA]
Justification:
There are no theoretical results included in our work.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have tried our best to include all the relevant details about every single experiment presented in this paper. We present the details of our experiments in SM Sections E and G These settings hold true for all experiments; unless otherwise stated.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
58

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: The codes and instructions on using them are available with our submission. The datasets are published in an anonymous repository and a URL. is provided. The codes will become public if the paper gets accepted.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ ptsips.cc/pulli//guidde Saape///ooeSuumiieeeeeee seriissssonpoiioio
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Please refer to SM Section Ablations of the architecture hyperparameters and the reasons for choosing them are provided in SM Sections F2 and
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: The most relevant statistical variability for this work is the randomness training. Due to the high computational costs of our trainings, we did not provide this information for the main experiments. However, we provide sensitivity to training randomness for a few datasets in SM Section E. We believe that the sensitivity generalizes to the rest of the experiments in work.
Guidelines:
The answer NA means that the paper does not include experiments.
59

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [No]
Justification: While we have not reported the execution time of individual experiments, we have provided the type of hardware and an estimate of total GPU hours for reproducing all the experiments SM Section
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make it into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http:///earrrsscccpppull/ccllcca
Answer: [Yes]
Justification: The main goal of the current work is to enhance ML models that surrogate computationally expensive physics simulators. Any misuse of the current work towards building harmful tools strongly condemned by the authors.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
60

---

Answer: [NA]
Justification: It is very hard to draw direct connection from the scope of this work to any societal impact.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (c.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: We fail cannot connect such misuses the physics applications of our work.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
The use of all existing code is properly credited. The source of all the existing datasets are mentioned and their licenses respected.
Guidelines:
The answer NA means that the paper does not existing assets.
61

---

The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re- re-packaged, both the original license and the license of the derived asset (if has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provide full instructions on how to use the data/code alongside each. The details of the datasets are fully explained the provided S.M.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve any form of crowdsourcing nor human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This work does not involve any form of crowdsourcing nor human subjects.
62

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question:
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The of LLMs in this work is only limited writing and editing.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our L.LM for what should or should not be described.
(https://neura) http://eerrrreeecco/o///e/e/////
63


Reviews:
[{"review_id": "4XSI0LEw5u", "reviewer": "Reviewer_1KFB", "summary": "This paper proposes  a novel graph-based framework for learning solution operators of PDEs on arbitrary, non-Cartesian domains. RIGNO uses an end-to-endGNN with a hierarchical \"regional mesh\" to achieve space-time resolution invariance and handle irregular point clouds. The framework introduces key innovations including regional mesh construction, edge masking, and temporal fractional pairing for robust and generalizable operator learning. Experiments on 13 unstructured and 8 Cartesian-grid PDE datasets show that RIGNO significantly outperforms prior methods (e.g., MeshGraphNet, FNO, CNO) in both accuracy and generalization.", "questions": "None.", "limitations": "None.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n1. The paper is novel to the best of my knowledge. It introduces a unique combination of graph neural operators and multi-scale regional meshes to handle non-Cartesian geometries with space-time resolution invariance, which is not sufficiently addressed in prior work.\n\n2. Writing is clear and well-structured.\n\n**Weakness**\n1. The current implementation focuses only on 2D domain. No empirical evidence or 3D experiments are shown.\n2. Fig 2 is small and not clear.\n3. Table1 looks strange. I believe the authors could find more baselines to fill the table.", "was_revised": false, "final_justification": "", "rebuttal": "We begin by thanking the reviewer for their thoughtful comments and constructive suggestions. We address their concerns in detail below.\n\n1. **Extension to 3D:**\nAs the reviewer rightly points out, the current work does not include empirical results on 3D problems. However, we would like to emphasize that all components of the proposed framework are readily extendable to 3D settings, and we do not anticipate any significant technical hurdles or computational limitations in doing so. In the graph construction stage (see Section \"Underlying Graphs\" on page 3 and SM Section A), all key concepts naturally generalize to 3D point clouds. The encoder and decoder edges defined in 2D using a support radius can be similarly defined in 3D, where the support sub-regions become spherical instead of circular. The processor edges and the support radii of the regional nodes can be constructed using tetrahedral meshes in 3D in place of 2D triangular meshes. Importantly, once the 3D input graphs are constructed, the underlying GNN architecture remains unchanged, requiring no architectural modifications to operate on 3D data.\n\n2. **Scalability:**\nOn the other hand, the 2D datasets used in our study are large-scale, with most containing 1K trajectories and over 28K input-output pairs, each defined on 16K spatial points. This level of scale significantly exceeds that of typical 2D datasets in the literature and is comparable to some 3D datasets, which typically include fewer than 1K input-output pairs at spatial resolutions ranging from 10K to 40K points. We believe that RIGNO’s strong performance on these challenging 2D benchmarks demonstrates both its scalability and its potential for extension to 3D problems. To further support this claim, we highlight the computational advantages of RIGNO over MeshGraphNet [1]. Notably, GNNs have already been employed successfully in 3D problems using MeshGraphNet and have been scaled to even larger problems as shown in [2]. To quantify this, we benchmarked the computational cost of RIGNO-18 against MeshGraphNet with 18 message-passing blocks on the Poisson-Gauss dataset (see table below). While MeshGraphNet saturates the memory of an NVIDIA GeForce RTX 3090 (24GB) with a batch size of 8, RIGNO (with 2.3K regional nodes) requires only half the memory, supporting a batch size of 16. When using fewer regional nodes (~0.5K), the batch size can be further increased to 32. Moreover, training and inference times are reduced to one-third compared to MeshGraphNet. These computational advantages do not come at the cost of performance. On the contrary, they lead to significant improvements, as evidenced by the results in Table 1 of the paper and SM Table F.2. Across many datasets, RIGNO achieves up to an order of magnitude reduction in test error compared to MeshGraphNet with similar configurations. Given that MeshGraphNet has been successfully tested on 3D problems, RIGNO's reduced computational/memory costs further support our claim that RIGNO can readily be extended to 3D problems. We aim to do so in future work. \n\n---\n> **Table 1: Computational cost of RIGNO and MeshGraphNet (MGN) with 18 message passing steps on the unstructured Poisson-Gauss dataset with 9216 input points. All experiments are done on a single NVIDIA GeForce 3090 (24GB). Training times are measured with 1024 training samples using the maximum batch size, and inference time corresponds to a single sample. RIGNO+ corresponds to a subsample factor of 16.0, resulting in only 576 regional nodes.**\n>\n> |  | Maximum batch size | Training time per epoch [s] | Inference time [ms] |\n> |:---:|---:|---:|---:|\n> | MGN | 8  | 53.1 | 25.5 |\n> | RIGNO | 16 | 28.4 | 12.7 |\n> | RIGNO+  | 32 | 17.9 | 7.2  |\n\n---\n\n\n6. **Further performance benchmarks:**\nWe benchmarked RIGNO and two baselines in terms of memory usage, training time, and inference time across various model configurations. The results indicate that RIGNO is consistently more memory-efficient than the baselines, while yielding higher accuracies and comparable inference times. We kindly refer the reviewer to point 6 in our response to Reviewer Z5X2 for a broader discussion regarding training and inference times, and point 6 in our response to Reviewer 9Nmc regarding memory usage. We will include these results, along with additional performance benchmarks, in a camera-ready version (CRV) if the paper is accepted.\n\n3. **Formatting issues:**\nWe thank the reviewer for their constructive suggestions regarding Figure 2 and Table 1, and will immediately correct these issues for the CRV, if accepted. Regarding Table 1, since the main advantage of the proposed architecture is for arbitrary domains, we have included many baselines that work on point clouds and only a few that work on uniform Cartesian grids. The extra page in a CRV allows us to simply separate the tables, especially since we plan to add Transolver as a baseline following the suggestion of reviewer AxpQ.\n\nWe hope that our responses have satisfactorily addressed all of the reviewer’s comments and concerns, and we kindly ask the reviewer to consider updating their assessment in light of these clarifications.\n\n*[1] Tobias Pfaff et al., Learning mesh-based simulation with graph networks. NeurIPS (2020).*\n\n*[2] Remi Lam et al., Learning skillful medium-range global weather forecasting. Science 382, 1416-1421 (2023).*"}, {"review_id": "TivEerF0EF", "reviewer": "Reviewer_Z5X2", "summary": "The authors propose RIGNO (Region Interacting Graph Neural Operator), an end-to-end graph neural network for training the solution operators of PDEs on arbitrary domains. The main motivation for the article is that many existing Neural Operators (NOs) assume that the underlying domain of the PDEs is discretized on a regular grid. The results demonstrate that RIGNO outperforms the baselines on both unstructured and structured grids.", "questions": "See the weaknesses", "limitations": "The authors discussed the method limitations.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n* RIGNO demonstrate better performance compare with state-of-art baselines across a wide and challenging range of benchmarks.\n* The paper introduces several novel techniques: \n\t* **edge masking**: the core idea is to randomly disable, or \"mask,\" a fraction of the edges in the graph before each message-passing step. This process is dynamic, with a different random set of edges being masked at every step and epoch.\n\t* **temporal fractional pairing** is a novel fine-tuning strategy designed to improve the generalization of time-continuous neural operators to time steps smaller than those seen during training.\n* The authors are also write about the limitations of their approach.\n* The introduction and related work sections clearly formulated the gap in the literature (the need for a multi-scale, resolution-invariant GNN operator) and position RIGNO as the solution.\n\n\n**Weaknesses**\n* The paper is purely empirical. It lacks theoretical justification for its key components:\n\t* Are there theoretical guarantees on why edge masking improve resolution invariance? It is presented only as an empirical result.\n* RIGNO contains a lot of hyperparameters. The article does not provide a sensitivity analysis for these hyperparameters. This complexity could make the model difficult to tune and apply for new problems, potentially hindering its practical application.\n* For a method to be a practical surrogate solver, its speed is as important as its accuracy. Therefore, it would be beneficial to include information in the article about the training and inference times for both RIGNO and the baselines.", "was_revised": true, "final_justification": "I get answer on my question. It helped me to better understand the paper.", "rebuttal": "We begin by thanking the reviewer for their thoughtful comments and constructive suggestions. We address their concerns in detail below.\n\n1. **Justification for edge masking:**\nAlthough we have not provided a formal theoretical analysis of edge masking in our manuscript, we have a clear and intuitive idea of how it contributes to resolution invariance. Considering the graph structure in the encoder, a regional node receives information from multiple physical nodes. With high-resolution inputs, many of these nodes are *redundant* in the sense that the underlying function can often be accurately reconstructed (e.g., via interpolated) using fewer nodes. During training, the model may implicitly ignore some of these redundant nodes and focus on a subset of dominant or central connections. At test time with different resolutions, these dominant nodes may no longer be present, leading to degraded performance. We identify this phenomenon as a key contributor to the resolution sensitivity observed in many GNN-based operator learning models and address it with the proposed method. By randomly masking edges during training, the model cannot consistently rely on the same high-correlation connections. In the absence of these edges, the model is forced to use the weaker correlations. This principle extends naturally to both the processor and decoder stages. Additionally, learning from weaker or secondary correlations in the processor enables the model to better integrate information across multiple spatial scales, rather than overfitting to dominant patterns in the training resolution.\n\n2. **Justification for fractional pairing:**\nSimilarly, the effectiveness of the fractional pairing strategy can be understood through a simple observation. In conventional training setups (including all2all pairing), the model is only exposed to lead times that are integer multiples of the dataset’s time resolution Dt. As a result, any smaller lead time (e.g., Dt/2) is out-of-distribution (OOD) for the model. During inference, this distribution shift can lead to noticeable performance degradation when the model is queried on shorter, unseen lead times. Our proposed fractional pairing fine-tuning strategy directly addresses this issue by explicitly extending the distribution of lead times seen during training. In simple terms, it increases the richness of our data distribution by approximating continuous time dynamics, thereby reducing generalization gaps and improving accuracy during inference.\n\n3. **Theoretical analysis:**\nAlthough mathematical formalization is lacking for much of the neural operator community, it is more straightforward to formalize the intuition behind edge masking. One approach to construct a formal definition could view edge masking as a Bayesian regularization on the input/output space of a GNN layer and follow a similar strategy as for dropout [1]. We thank the reviewer for this thoughtful suggestion and recognize it as a valuable direction for future research. We also appreciate the reviewer’s observation regarding the lack of detailed justification of the effectiveness of these methods, and we will gladly include a more comprehensive discussion in the camera-ready version (CRV), if the paper is accepted. However, providing rigorous theoretical explanations or guarantees would require significant additional effort and is beyond the scope of the current work.\n\n4. **Many hyperparameters:**\nApart from standard optimization hyperparameters common to most machine learning frameworks, the size and performance of RIGNO are primarily sensitive to just **two key architectural hyperparameters**: the number of processor blocks and the latent dimension of edge/node features. In fact, one of the core design principles of our framework is to minimize the need for extensive hyperparameter tuning by using **a fixed set of hyperparameters across all tasks**. To this end, we identified a robust configuration by tuning on a single internal dataset and applied the same hyperparameters across all other datasets without further adjustment. As shown in Table 1 of the submission, both the proposed RIGNO-18 and RIGNO-12 architectures deliver consistently strong performance across diverse benchmarks, in contrast to many of the baselines that require task-specific tuning. The complete list of hyperparameters is provided in SM Table E.1, and this configuration can be readily applied to new problems or datasets. Note that the only difference between RIGNO-18 and RIGNO-12 lies in the number of processor blocks. While we report all hyperparameters for transparency, most of them do not require tuning. If the paper is accepted, we will clearly emphasize this point in the CRV.\n\n5. **Sensitivity to hyperparameters:**\nWe have analyzed the sensitivity of RIGNO with respect to both accuracy and inference time across four hyperparameters in SM Section F.3 (see SM Figure F.1). In particular, we demonstrate that overlap factors and the number of regional mesh levels do not require tuning. Based on our analysis, we recommend using minimal overlaps, as they offer limited accuracy gains relative to their computational cost, and maximal regional mesh levels, which provide substantial accuracy improvements with minimal added cost. The number of regional nodes is another important hyperparameter, and is studied in SM Section F.2 (see SM Table F.2). Our results show that while using a coarser regional mesh can improve efficiency, it reduces model capacity. We find that using approximately 4000 regional nodes provides a good trade-off between accuracy and efficiency. While the optimal number may vary depending on the target function, this configuration has consistently yielded strong results across all datasets.\n\n6. **Performance benchmarks:**\nBelow, we provide two tables reporting the training and inference times for GINO, RIGNO, and MeshGraphNet (MGN) [1] across various model configurations, all evaluated on a single NVIDIA GeForce RTX 3090 GPU (24GB). GINO is included due to its consistently accurate performance across benchmarks, while MGN is selected for its architectural similarity to RIGNO. Owing to the limited time available for evaluation and the consistent accuracy advantages of RIGNO over MGN (see Table 1 in the main text), we report test errors only for RIGNO and GINO in the tables below. For training, we use the maximum possible batch size for RIGNO and MGN, and a batch size of 64 for GINO. In scenarios involving varying geometries across training samples, such as in the Elasticity experiment, the original implementation of GINO does not support batch processing, requiring us to use a reduced batch size of 1. We employ 18 message-passing blocks in the processor for both RIGNO and MGN, and evaluate configurations with 96 and 128 latent features. For GINO, we report results for two variants consisting of 3 and 7 FNO layers, respectively. The results show that the training times of RIGNO are slower than those of GINO for datasets where the geometry and input coordinates are fixed between training samples, and are faster for datasets with variable geometries, mainly due to the support of batch processing. We acknowledge the former case is a disadvantage of the proposed model; however, as training time is a fixed initial cost, the difference in speed is more critical at inference time. In this setting, RIGNO and GINO have very similar performance, yet RIGNO is still able to significantly surpass GINO's performance in terms of accuracy. While this subset of experiments is necessarily limited due to the time constraints of the rebuttal period, we believe the reported metrics are representative and consistent across settings. Admittedly, longer training times are a drawback; however, we view them as a worthwhile tradeoff for the consistent gains in accuracy and robustness achieved across a diverse range of physical systems with RIGNO over GINO (Main Text Table 1). Importantly, inference times, which are far more relevant for practical deployment, remain comparable across all models considered in this rebuttal and in the main text (Table 1). Moreover, our benchmarks show that RIGNO is up to 3x more memory-efficient than GINO. Please see point 6 in our response to Reviewer 9Nmc. Thanks to its multi-scale architecture and computational efficiency, RIGNO emerges as a powerful end-to-end GNN that achieves state-of-the-art accuracy and robustness, while remaining competitive with strong non-GNN baselines in terms of speed and memory usage. We believe this work reinforces the viability of GNNs for solving PDEs and opens the door to future developments aimed at further reducing training costs. We thank the reviewer for raising this point and will present these results along with complementary performance benchmarks in a CRV, if accepted.\n\n---\n>**Table 1: Performance benchmarks on the unstructured Poisson-Gauss. RIGNO+ corresponds a subsample factor of 16.0, resulting in only 576 regional nodes.**\n>\n>||Parameters|Training time per epoch [s]|Test error [%]|Inference time [ms]|\n>|:---|:---:|:---:|:---:|:---:|\n>|GINO (3)|3.6|**1.33**|5.24|7.8|\n>|GINO (7)|8.4|**1.72**|8.01|10.1|\n>|MGN (96)|1.4|43.3|n/a|16.1|\n>|MGN (128)|2.5|53.1|n/a|25.5|\n>|RIGNO (96)|1.4|24.9|**2.23**|8.1|\n>|RIGNO (128)|2.5|28.4|**2.09**|12.7|\n>|RIGNO+ (96)|1.4|15.5|4.19|**5.2**|\n>|RIGNO+ (128)|2.5|17.9|3.25|**7.2**|\n\n>**Table 2: Performance benchmarks on the Elasticity dataset.**\n>\n>||Parameters|Training time per epoch [s]|Inference time [ms]|\n>|:---|:---:|:---:|:---:|\n>|GINO (3)|3.6|18.2|4.1|\n>|GINO (7) |8.4|27.1|6.4|\n>|RIGNO (96)|1.4|**2.4**|**1.5**|\n>|RIGNO (128)|2.5|**3.1**|**1.9**|\n---\n\nWe hope that our responses have satisfactorily addressed all of the reviewer’s comments and concerns, and we kindly ask the reviewer to consider updating their assessment in light of these clarifications.\n\n*[1] Pierre Baldi and Peter Sadowski, Understanding dropout. NeurIPS (2013).*"}, {"review_id": "DddBRJhl25", "reviewer": "Reviewer_9Nmc", "summary": "This paper presents the Region Interacting Graph Neural Operator (RIGNO), a multi-scale, resolution-invariant graph neural operator that builds upon established architectural principles while introducing several interesting enhancements. The work follows a solid foundation using latent space representations and adopts the encoder-processor-decoder structure familiar from MeshGraphNet, enhanced with a novel multiscale regional mapping component.", "questions": "Based on my understanding, the paper proposes to use graph neural operators for the flexibility of graph. Which component from RIGNO makes it signficiantly outperform previous state-of-the-art model such as CNO and Poseidon/SCOT on regular grids? Are these components applicable to other models?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strengths\nThe paper demonstrates several commendable qualities. The experimental evaluation is particularly thorough, encompassing a comprehensive suite of PDE test problems and comparing against a diverse set of baseline methods. The reported performance improvements are impressive and suggest the method's practical value. The incorporation of edge masking for regularization and temporal processing through adaptive layer normalization with gating mechanisms shows thoughtful engineering.\n\n## Weaknesses and Questions\nWhile the work shows strong empirical results, I believe it could benefit from deeper analysis in a few key areas:\n- Novelty Clarification: The paper skillfully combines several existing techniques, but it would strengthen the contribution to more clearly articulate which specific innovations drive the state-of-the-art performance. A more explicit discussion of the novel aspects beyond the multiscale regional mapping would help readers understand the core contributions.\n- Ablation study: The work would significantly benefit from ablation studies examining each major component's contribution to overall performance. This would help identify which design choices are most critical and provide valuable insights for future work in this area.\n- Computational Efficiency: Including runtime and memory usage comparisons would provide important practical insights, especially given the multi-scale nature of the approach, which readers would find valuable for implementation considerations.", "was_revised": false, "final_justification": "", "rebuttal": "We begin by thanking the reviewer for their thoughtful comments, constructive suggestions, and their recognition of the value of our work. We address their concerns in detail below.\n\n1. **Contributing factors to RIGNO's performance:**\nAmong the baselines presented in Table 1, most (excluding CNO and MeshGraphNet [1]) rely on global interactions; either through convolutions in the Fourier space or multi-head attention mechanisms. In contrast, we attribute the state-of-the-art accuracy and robustness of RIGNO primarily to the high expressive power of *local* message passing for fine-grained manipulation and propagation of information. This advantage is especially prominent in time-dependent and advection-dominated problems, whereas for elliptic problems (AF, Elasticity, and Poisson-Gauss), where global interactions are more crucial, the gains are less pronounced. Compared to MeshGraphNet [1], which also utilizes message passing, RIGNO introduces several key architectural components that enhance both efficiency and expressivity. The introduction of regional nodes facilitates multi-scale aggregation of input data and yields significant computational benefits (see SM Table F.2 and Point 2 in our response to Reviewer 1KFB). While we benchmarked MeshGraphNet with the same number of message-passing blocks, SM Figure F.1 indicates that such models benefit substantially from deeper message-passing pipelines. RIGNO enables this depth efficiently by performing message passing on coarse regional graphs, which are less costly to process. To mitigate the limitations of purely local interactions, RIGNO incorporates *multi-scale edges* in the processor. These long-range connections significantly enhance accuracy, without a notable increase in computational cost, as shown in the right-most plot of SM Figure F.1. We appreciate the reviewer’s suggestion to provide a broader discussion of the novel components driving RIGNO’s performance, and will include this analysis in the camera-ready version (CRV), if accepted. Below, we further elaborate on two of our key contributions, *edge masking* and *fractional pairing*, whose design and benefits will also be described in greater detail in the CRV. We emphasize that both techniques are framework-agnostic and easily transferable. Edge masking can be integrated into any GNN-based PDE solver to improve generalization across spatial resolutions. Similarly, fractional pairing provides a lightweight, effective fine-tuning strategy for mitigating temporal generalization gaps and can easily be **employed in any operator learning framework**.\n\n2. **Intuitions behind edge masking:**\nConsidering the graph structure in the encoder, a regional node receives information from multiple physical nodes. With high-resolution inputs, many of these nodes are *redundant* in the sense that the underlying function can often be accurately reconstructed (e.g., via interpolation) using fewer nodes. During training, the model may implicitly ignore some of these redundant nodes and focus on a subset of dominant or central connections. At test time with different resolutions, these dominant nodes may no longer be present, leading to degraded performance. We identify this phenomenon as a key contributor to the resolution sensitivity observed in many GNN-based operator learning models and address it with the proposed method. By randomly masking edges during training, the model cannot consistently rely on the same high-correlation connections. In the absence of these edges, the model is forced to use the weaker correlations. This principle extends naturally to both the processor and decoder stages. Additionally, learning from weaker or secondary correlations in the processor enables the model to better integrate information across multiple spatial scales, rather than overfitting to dominant patterns in the training resolution.\n\n3. **Intuitions behind fractional pairing:**\nSimilarly, the effectiveness of the fractional pairing strategy can be understood through a simple observation. In conventional training setups (including all2all pairing), the model is only exposed to lead times that are integer multiples of the dataset’s time resolution Dt. As a result, any smaller lead time (e.g., Dt/2) is out-of-distribution (OOD) for the model. During inference, this distribution shift can lead to noticeable performance degradation when the model is queried on shorter, unseen lead times. Our proposed fractional pairing fine-tuning strategy directly addresses this issue by explicitly extending the distribution of lead times seen during training. In simple terms, it increases the richness of our data distribution by approximating continuous time dynamics, thereby reducing generalization gaps and improving accuracy during inference.\n\n4. **Ablation studies:**\nWe have conducted a comprehensive ablation study to evaluate the core components of the proposed framework and to identify their individual contributions to overall performance. While we present the most critical results in the main text, additional ablations and supporting evidence are provided in the Supplementary Materials (SM). In Figure 3a and the accompanying discussion on page 7, we ablate the proposed *edge masking* strategy, highlighting its role in enhancing resolution invariance. We also briefly discuss its computational advantages. Further, in SM Section F.4, we demonstrate the broader applicability of edge masking by incorporating it into MeshGraphNet [1], where it similarly leads to a substantial reduction in the test error with unseen resolutions. In both figures, it is evident that edge masking also reduces the test error with the training resolution by facilitating more effective optimization. The impact of the fractional pairing strategy and various time-marching schemes is ablated in Figures 3b and 3c. We also ablate the regional aggregation mechanism in the encoder and decoder in SM Section F.2 (see SM Table F.2). These results confirm that regional aggregation is a key contributor to RIGNO's performance. Notably, the trend in decreasing test error breaks when regional aggregation is omitted (at 16,384 regional nodes, same as the number of input points), underscoring its necessity for accurate learning. Additionally, in SM Section F.3 (see SM Figure F.1), we analyze the sensitivity of RIGNO’s accuracy and inference time to four major hyperparameters. From this analysis, we observe that neither the overlap factors nor the number of regional mesh levels require tuning: minimal overlaps are recommended due to their limited benefit relative to their computational cost, while maximal mesh levels are encouraged given their accuracy gains with negligible overhead. Moreover, we find that the most effective way to scale RIGNO’s capacity is through increasing either the number of processor blocks or the latent dimensionality of node and edge features. Both of these options yield consistent improvements in performance. We believe that all essential aspects of the framework have been carefully studied and reported. Other hyperparameters listed in SM Table E.1 (e.g., the number of hidden layers in MLPs) are included for completeness, though we find they have minimal impact on model accuracy or generalization.\n\n5. **Training/inference time benchmarks:**\nWe benchmarked GINO, RIGNO, and MeshGraphNet [1] in terms of training and inference times across various model configurations. The results indicate that RIGNO’s training times are consistently slower than those of GINO on datasets with fixed geometry and consistent input coordinates across training samples. We acknowledge this as a limitation of our model; however, since training is a one-time cost, inference speed is ultimately more critical for deployment. In this regard, RIGNO and GINO exhibit comparable inference times, while RIGNO consistently achieves higher accuracy. Moreover, unlike GINO, RIGNO supports batch processing for datasets with variable geometries (e.g., Elasticity), where it is also faster at training. Due to space constraints, we kindly refer the reviewer to point 6 in our response to Reviewer Z5X2 for a broader discussion. We thank the reviewer for highlighting this point and will include these results, along with additional performance benchmarks, in a CRV if the paper is accepted.\n\n6. **Memory usage benchmarks:**\nFollowing the reviewer’s suggestion, we evaluated the memory usage of RIGNO and the best-performing baseline, GINO, in terms of peak inference memory usage on the Elasticity and Poisson-Gauss datasets using a standard 24GB NVIDIA GPU. The results are summarized in the table below. For the Poisson-Gauss dataset, GINO requires approximately 370MB per sample, which limits the maximum batch size to 66. In contrast, RIGNO uses only 120MB per sample, allowing batch sizes up to 200. A similar trend is observed on the Elasticity dataset, where the smaller problem size permits even larger batches for both models. If the paper is accepted, we will gladly include these results, along with a comparison to additional baselines, in a CRV.\n\n---\n>**Table 1: Peak memory usage of RIGNO and GINO on the unstructured Poisson-Gauss and Elasticity datasets with a NVIDIA GeForce RTX 3090 (24GB). The two best results are highlighted at each column.**\n>\n>||Parameters [M]|Poisson-Gauss|Elasticity|\n>|:---|:---:|:---:|:---:|\n>|GINO (3)|3.6|370.0MB|80.3MB|\n>|GINO (7) |8.4|406.1MB|116.4MB|\n>|RIGNO (96)|1.4|**120.9MB**|**50.6MB**|\n>|RIGNO (128)|2.5|**159.3MB**|**64.6MB**|\n---\n\nWe hope to have addressed all of the reviewer’s comments and concerns in this response, and we kindly ask the reviewer to consider updating their assessment in light of these clarifications.\n\n*[1] Tobias Pfaff et al., Learning mesh-based simulation with graph networks. NeurIPS (2020).*"}, {"review_id": "2xE067pagS", "reviewer": "Reviewer_AxpQ", "summary": "This paper presents RIGNO, a new graph-based neural operator framework for solving PDEs on arbitrary domains. RIGNO proposes a new approach for encoding sets of points via a regional mesh approach, which enable to compress the physical information. It then uses a message passing architecture to predict future time-steps via an all2all training approach. This allows to generalize to different timestep intervals. To enable continuous timestep evaluation, they propose temporal fractional pairing to allow generalization to intermediate timesteps not accessible in the training set. Their approach is evaluated on a range of PDEs and exhibit competitive or sota performance.", "questions": "Did you try to test your model on out-of-domain data, e.g. on new PDE coefficients? As you give the pde coefficients as input, it could be interesting to check how the model behaves on unseen coefficients.\n\nIn the all2all training, you only provide future timesteps to predict right? Could the model actually learns to solve an inverse problem ? e.g. predicting the past frames?", "limitations": "Limitations have been mentioned. Negative societal impacts have not been mentioned in the work.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\nI think the paper is clearly written, in particular, there are very nice illustrations in appendix that clearly help for the comprehension of the paper. \n\nThe method has been evaluated on a large range of PDEs with substantial experiments. On all datasets, Rigno-18 outperforms the considered baselines. the proposed framework is a simple but effective encode-process-decode approach, which allows to process arbitrary geometries. \n\nI think the paper also tackles an important problem: generalization on arbitrary geometries. The authors could other important works that differ from graphs or transformers such as INRs, which is another direction for learning on arbitrary geometries [1, 2].\n\nI also appreciated the ablations studies and additional experiments provided, notably for the time processing component. \n\nWeaknesses:\n\nRegarding the novelty of the paper, the method is relatively simple and follows a similar direction than existing approaches. It follows the encode-process-decode approach that is particularly effective for graph-based neural operators and make use of the training strategy proposed in Poseidon. Their contributions include the edge masking and the Fractional pairing strategy which are relatively simple. However, they provide ablations showing the advantage of their incorporation.\n\nFor the ablations on the dataset and model size, it would have been nice to test until when model size or data size reaches a plateau. It has been on relatively small/medium setups only.\n\nThere are other important baselines that could have been used such as [1] and [3] for instance.\n\nAs you mentioned in the introduction, efficiency is a key property for learning efficient neural operators on arbitrary geometries. Graphs are known to be particularly costly. It is important to add additional experiments to show how RIGNO behaves in terms of computation with respect to considered baselines.\n\n[1] Serrano et al. Operator Learning with Neural Fields: Tackling PDEs on General Geometries, NeurIPS 2023.\n\n[2] Wang et al. GridMix: Exploring Spatial Modulation for Neural Fields in PDE Modeling, ICLR 2025.\n\n[3] Wu et al. Transolver: A Fast Transformer Solver for PDEs on General Geometries, ICML 2024.", "was_revised": true, "final_justification": "Most of my concerns have been addressed by the authors. Thus, I have upgraded my score to a borderline accept. I think the novelty of the paper is limited, but the provided experiments are consequent.\n\n- the authors provided the computational efficiency and acknowledge it as a limitation of the paper.\n- the INR literature will be discussed and CORAL has been compared to RIGNO.\n- they did out-horizon experiments for OOD.\n\nThey proposed to do OOD evaluation on unseen PDE parameters. I think it could be valuable.", "rebuttal": "We would like to thank the reviewer for their excellent and constructive suggestions. We hope to address all their concerns in detail below.\n\n1. **Novelty and impact**:\nWe appreciate the reviewer’s observation that the proposed contributions, edge masking and fractional pairing, are grounded in simple and intuitive ideas. Nonetheless, these play a critical role in RIGNO's ability to generalize to unseen spatial and temporal resolutions. We see the relative simplicity of these approaches as an advantage, as these methods are easily transferable across frameworks. Thus, this work contributes a broader impact in the operator learning community than a standalone neural operator. For example, edge masking can be readily integrated into any GNN-based PDE solver, while fractional pairing offers a lightweight and practical strategy for addressing generalization across temporal resolutions in **any** operator learning framework.\n\n2. **Dataset and model scaling experiments:**\nPerformance plateaus can arise from either a large generalization gap (i.e., insufficient data) or limited model expressivity (i.e., small model capacity). In Figure 2, our primary aim is to show that both issues can be effectively mitigated within the proposed framework; by increasing the amount of data in the former case, and scaling the model size in the latter. Specifically, we used 1024 training samples for the time-independent datasets and 1024 trajectories for the time-dependent ones (7168 input-output pairs). For single-task operator learning, this represents a relatively data-rich regime, especially considering the high resolution of the input/output functions. In this setting, we do not expect to observe a plateau unless the model becomes exceedingly large. To further support this within the limited time of the rebuttal period, we conducted a small-scale experiment using only 128 samples/trajectories from two datasets. The results (see table below) indicate that although the training error continues to decrease (overfit), the test error plateaus around 1.8M parameters. Linearly extrapolating this result to the training conditions of the paper (outlined in Fig. 2), suggests that a plateau would occur at approximately 16M parameters or more. As the reviewer rightfully suggests, a more extensive version of these experiments would be a valuable addition. We will add these to a camera-ready version (CRV), if accepted.\n\n\n---\n>**Table 1: Scaling behavior of test and training errors (%) with model size for Poisson-Gauss (PG) and NS-Gauss (NS) datasets, using 128 training samples for PG and 128 trajectories for NS.**\n>\n>| Parameters [M] | PG Test | PG Train | NS Test | NS Train |\n>|:--------------:|--------:|---------:|--------:|---------:|\n>| 0.1            | 19.93   | 13.06    | 39.28   | 31.46    |\n>| 0.5            | 16.77   | 9.05     | 31.03   | 17.39    |\n>| 1.1            | 13.73   | 5.97     | 15.29   | 11.12    |\n>| 1.8            | 13.14   | 4.56     | 13.01   | 9.57     |\n>| 2.8            | 13.45   | 3.99     | 10.79   | 8.21     |\n>| 4.1            | 13.23   | 3.18     | 11.46   | 7.14     |\n\n---\n\n\n3. **Alternative approach:** INRs offer an alternative approach for learning on arbitrary geometries by encoding both the input and output functions into single latent feature vectors, with the mapping learned in this latent space. Given both the popularity and effectiveness of transformers and graph-based models, as well as RIGNO's position as an advancement among these methods, we have focused on these baselines. Nonetheless, we would like to point out that among the 8 models used for comparison, we have used the FNO and CNO, which are not graph or transformer based, as well as the GeoFNO and FNO DSE, two popular approaches for data on general geometries. These baselines cover graph, convolution, attention, and Fourier-based frameworks. Nonetheless, since INRs also address one of the core problems considered in RIGNO, we agree that they warrant discussion in the broader context of operator learning. We thank the reviewer for this suggestion and will include a brief discussion of INRs in the CRV, if the paper is accepted.\n\n4. **Further benchmarks (Transolver):**\nBased on the reviewer’s helpful suggestion, we evaluated Transolver as an additional baseline. However, we found that this model scales poorly, severely hindering its ability to learn on large-scale datasets. This limitation stems from its use of global interactions, implemented via independent MLPs acting between each input point and a set of sparse tokens. The resulting large number of MLP evaluations leads to substantial memory consumption, limiting the model size to at most 0.5M parameters on an RTX 4090 GPU. We trained Transolver on two representative datasets, CE-RP and Wave-C-Sines, and observed relatively high test errors of 17.8% and 27.7%, respectively. These are substantially higher than the corresponding RIGNO-12 errors of 4.92% and 6.25%, achieved under similar computational constraints. We thank the reviewer for this valuable suggestion and will include Transolver results for all datasets in the CRV, if the paper is accepted.\n\n5. **Comparison with CORAL:**\nFollowing the reviewer's suggestion, we investigated RIGNO's performance in comparison to CORAL. Given the time constraints and major differences in the training strategies, we do not train CORAL, and rather rely on the test results in Table 3 of [1]. The NACA-Euler dataset of [1] is the same airfoil dataset that we use in RIGNO (AF), with the difference that the target variable used there is the Mach number. We conducted new experiments with RIGNO with the settings of [1] (1000 training samples, 200 test samples), and used the same metric as in [1] (relative L2 error without normalization) to obtain the following results in Table 2, where we also report the test errors for the second-best and third-best baselines in [1]. Based on the results in Table 1 of the submission, RIGNO's accuracy is most striking for hyperbolic and parabolic problems, and does not have a large margin with the other baselines on elliptic problems. Nonetheless, RIGNO-18 achieves 0.35% error on 200 test samples, outperforming CORAL (0.59%) and the other baselines considered in [1] by a considerable margin.\n\n---\n>**Table 2: Relative L2 errors (\\%) without normalization on 200 test samples. All models are trained on 1000 samples.**\n>\n>|  | NACA-Euler |\n>|:--------------:|-----------:|\n>| RIGNO-18       | 0.35  |\n>| RIGNO-12       | 0.40    |\n>| CORAL          | 0.59   |\n>| Factorized-FNO | 0.62   |\n>| Geo-FNO        | 1.58       |\n---\n\n\n6. **Further performance benchmarks:**\nWe benchmarked GINO, RIGNO, and MeshGraphNet [1] in terms of training and inference times across various model configurations. The results indicate that RIGNO’s training times are consistently slower than those of GINO on datasets with fixed geometry and consistent input coordinates across training samples. We acknowledge this as a limitation of our model; however, since training is a one-time cost, inference speed is ultimately more critical for deployment. In this regard, RIGNO and GINO exhibit comparable inference times, while RIGNO consistently achieves higher accuracy. Moreover, unlike GINO, RIGNO supports batch processing for datasets with variable geometries (e.g., Elasticity), where it is also faster at training. Our benchmarks show that RIGNO is also consistently more memory-efficient than GINO, enabling it to be inferred with ~3x more samples. Due to space constraints, we kindly refer the reviewer to point 6 in our response to Reviewer Z5X2 for a broader discussion regarding training and inference times, and point 6 in our response to Reviewer 9Nmc regarding memory usage. We thank the reviewer for highlighting this point and will include these results, along with additional performance benchmarks, in a CRV if the paper is accepted.\n\n7. **Out-of-domain experiments:**\nWe explore out-of-domain performance in the form of solutions at **unseen time steps**. Figure 3c in the main text shows the performance of RIGNO on NS-SL for time steps 0-14, seen during training, as well as for *unseen time steps* 16-20. These later time steps have stronger vortices and different dynamics, yet the error does not diverge from the general trend, indicating strong out-of-domain performance. SM Section F.5 and SM Figure F.3 offer a broader evaluation and discussion regarding this test with the NS-PwC dataset. For variation in the parameters of the PDEs, there exist several datasets which include *out-of-distribution* test sets, available in [2]. Given the limited time of the rebuttal period, we have been unable to perform these experiments; however, we will gladly add these experiments to the CRV, if accepted.\n\n8. **Inverse problems:**\nIn all2all training, the reviewer is correct: we only provide future time steps. Neural operators can learn arbitrary maps between parameters; however, inverse maps are often ill-posed. For this reason, approaches for learning inverse problems typically use an invertible architecture, but this may impose challenges in learning expressive maps with high-frequency features [3]. We believe studying the use of neural operators for inverse problems deserves a thorough investigation beyond the scope of the current work, and we reserve this task for future work. \n\nWe hope that our responses have satisfactorily addressed all of the reviewer’s comments and concerns, and we will happily provide any additional clarifications. We kindly request the reviewer to consider updating their assessment in light of these remarks.\n\n*[1] Serrano et al, Operator Learning with Neural Fields: Tackling PDEs on General Geometries. NeurIPS 2023.*\n\n*[2] Raonic et al, Convolutional Neural Operators for robust and accurate learning of PDEs. NeurIPS 2023.*\n\n*[3] Long et al, Invertible Fourier Neural Operators for Tackling Both Forward and Inverse Problems. PMLR 2024.*"}]

Author rebuttal:

Meta-review: This paper presents RIGNO, a graph-based framework for robust and accurate operator learning of PDEs on arbitrary domains. The reviewers initially highlighted concerns regarding clarity of the methodology, insufficient discussion of efficiency and scaling, and the need for stronger baselines. While they found the core idea promising, they felt that these issues limited the paper’s impact in its original form.

The authors’ rebuttal directly addressed these points. They provided detailed clarifications on the role of fractional pairing and other design choices, added comparisons against additional baselines, and included new efficiency benchmarks that demonstrated competitive training and inference performance. Reviewers acknowledged that these responses resolved their questions, improved their understanding of the method, and in some cases led to higher confidence in their evaluations.

Overall, the paper now presents a clear, well-supported framework with strong experimental evidence and thorough justification of design choices. With the reviewers’ concerns adequately resolved in the rebuttal, I recommend acceptance.
Decision: Accept (poster)