Space and time continuous physics simulation from partial observations
Authors: Steeven JANNY
Madiha Nadri
Julie Digne
Christian Wolf
Venue: ICLR 2024 spotlight (2024)
Keywords: Physics
simulation
interpolation
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Modern techniques for physical simulations rely on numerical schemes and mesh-refinement methods to address trade-offs between precision and complexity, but these handcrafted solutions are tedious and require high computational power. Data-driven methods based on large-scale machine learning promise high adaptivity by integrating long-range dependencies more directly and efficiently. In this work, we focus on computational fluid dynamics and address the shortcomings of a large part of the literature, which are based on fixed support for computations and predictions in the form of regular or irregular grids. We propose a novel setup to perform predictions in a continuous spatial and temporal domain while being trained on sparse observations. We formulate the task as a double observation problem and propose a solution with two interlinked dynamical systems defined on, respectively, the sparse positions and the continuous domain, which allows to forecast and interpolate a solution from the initial condition. Our practical implementation involves recurrent GNNs and a spatio-temporal attention observer capable of interpolating the solution at arbitrary locations. Our model not only generalizes to new initial conditions (as standard auto-regressive models do) but also performs evaluation at arbitrary space and time locations. We evaluate on three standard datasets in fluid dynamics and compare to strong baselines, which are outperformed in classical settings and the extended new task requiring continuous predictions.

Published as a conference paper at ICLR 2024
SPACE AND TIME CONTINUOUS PHYSICS SIMULATION FROM PARTIAL OBSERVATIONS
Steeven Janny
Madiha Nadri
LIRIS, INSA Lyon, France steeven. janny nsa-lyon. fr
LAGEPP, Univ. Lyon 1, France madiha hawriaaaaaddooooooo féuniv-lyon! fr
Julie Digne
Christian Wolf
LIRIS, CNRS, France
Naver Labs Europe, France christian, www.enaverrlbber. com
julie. digne@cnrs.
ABSTRACT
Modern techniques for physical simulations rely on numerical schemes and mesh- refinement methods to address trade-offss between precision and complexity, but these handcrafted solutions are tedious and require high computational power. Data-driven methods based on large-scale machine learning promise high adap- tivity by integrating long-range dependencies more directly and efficiently. In this work, we focus on fluid dynamics and address the shortcomings of large part of the literature, which are based on fixed support for computations and predictions in the form of regular or irregular grids. We propose novel setup perform pre- dictions in a continuous spatial and temporal domain while being trained on sparse observations. We formulate the task as a double observation problem and propose a solution with two interlinked dynamical systems defined on, respectively. the sparse positions and the continuous domain, which allows to forecast and interpo- late a solution from the initial condition. Our practical implementation involves recurrent GNNs and spatio-temporal attention observer capable of interpolating the solution at arbitrary locations. Our model not only generalizes to new initial conditions (as standard rugo-segressive models do) but also performs evaluation at arbitrary space and time locations. We evaluate on three standard datasets in fluid dynamics and compare to strong baselines, which are outperformed both classical settings and in the extended new task requiring continuous predictions.
INTRODUCTION
The Lavoisier conservation principle states that changes in physical quantities in closed regions must be attributed to either input, output, or source terms. By applying this rule at an infinitesimal scale, retrieve partial differential equations (PDEs) governing the evolution of large majority of physics scenarios. Consequently, the development of efficient solvers is crucial in various domains involving physical phenomena. While conventional methods (c.g. finite difference or finite volume methods) showed early success in many situations, numerical schemes suffer from high computa- tional complexity, particular for growing requirements on fidelity and precision. Therefore, there is a need for faster and more versatile simulation tools that are reliable and efficient, and data-driven methods offer a promising opportunity.
Large-scale machine learning offers natural solution to this problem. In this paper, we address data-driven solvers for physics, but with additional requirements on the behavior of simulator:
RJ. Data-driven- the underlying physics equation is assumed to be completely unknown. This includes the PDE, but also the boundary conditions. The dynamics must be discovered from finite dataset of trajectories, i.e. a collection of observed behaviors from the physical system,
R2. Generalization the method must be capable of handling new initial conditions that do not explicitly belong to the training set, without re-training fine-tuning.

---

Published as a conference paper at ICLR 2024
R3. Time and space continuous the domain of the predicted solution must be continuous in space and timd so that it can be queried at any arbitrary location within the domain of definition.
These requirements common in the field but rarely addressed altogether. R1 allows for handling complex phenomena where the exact equation might be unknown, and R2 supports the growing need for faster simulators, which consequently must handle new ICs. Space and time continuity (R3) are also useful properties for standard simulations since the solution can be made as fine as needed in certain complex areas.
This task requires learning from sparsely distributed observations only, and without any prior knowl- edge on the PDE form. In these settings, a standard approach consists of approximating the behav- ior of a discrete solver, enabling forecasting in an auto-regressive fashion Pfaff et al. (2020); (Jann):Jannny al. (2023); Sanchez-Gonzzlez et al.] (2020) losing therefore spatial and temporal continuity. In- deed, auto-regressive models assume strong regularities in the data, such as static spatial lattice and uniform time steps. For these reasons, generalization to new spatial locations or intermediate time steps is not straightforward. These methods satisfy RI and R2, but not R3. In another trend, Physics-Infornee Neural Networks (PINNs) learn a solution on a continuous domain. They leverage the PDE operator to optimize the weights of a neural network representing the solution, and cannot generalize to new ICs, thus violating RI and R2.
In this paper, we address RI, R2 and R3 altogether in new setup involving two joint dynamical systems. RI and R2 are satisfied using an auto-regressive discrete-time dynamics learned from the sparse observations and producing trajectory in latent space. Then, R3 is achieved with state observer derived from a second dynamical system in continuous time. This state observer relies on transformer-based cross-attention to enable evaluation arbitrary spatio-temporal locations. In nutshell: (a) We propose a new setup to address continuous space and time simulations of physical systems from sparse observation, leveraging insights from control theory. (b) We provide strong theoretical results indicating that our setup is well-suited to address this task compared to existing baselines, which are confirmed experimentally on challenging benchmarks. (c) We provide exper- imental evidence that our state observer is more powerful than handcrafted interpolations for the targeted task. (d) With experiments on three challenging standard datasets (Navier Yin et al. (2022);
Stokes (2009), Shallow Water Yin et al. (2022); Galewsky et al.] (2004), Eagle Janny et al. (023)) and against state-of-the-arr methods (MeshGraphNet [MGf) et al. (2020). DINO Yin et al. (022),, MAgNet (Boussif et al.) (2022)), we show that our results generalize to wider class of problems, with excellent performances.
2 RELATED WORKS
Autoregressive models have been extensively used to replicate the behavior of iterative solvers in discrete time, especially in cases where the PDE is unknown or generalization new initial condi- tions is needed. These models come in various internal architectures, including convolution-basee models for systems observed on a dense uniform grid (Stachenfeld et al. 2021 Guen & Thome 2020 Bézenac et al. 2019) and graph neural networks (Battaglia et al.] (016)) that can adapt to ar- bitrary spatial discretizations Sanchez-Gonzalee
et al. 2020 Janny et al. [202a Li et al. 2018).
Such models have demonstrated a remarkable capacity to produce highly accurate predictions and generalize over long prediction horizons, making them particularly suitable for addressing complex problems such as fluid simulation (Piaff et al. [2000] Hann et al. 2021 Janny et al. (023). However, auto-regressive models are inherently limited to a fixed and constant spatio-temporal discretization grid, hindering their capability to evaluate the solution anywhere and at any time. Neural ordinary differential equations (Neural ODE Chen et al. (2018); : Dupont et al. (2019) offer a countermea- sure to the fixed timestep constraint by learning continuous ODEs discrete data using explicit solver, such as Fuler or Runge-Kurt Kutta methods. In theory. this enables the solution to be evaluated at any temporal location but in practice still relies on the discretization of the time variable. Moreover, extending this approach to PDEs is not straightforward. Contrarily to these approaches, we leverage the auto-regressive capacity and accuracy while allowing arbitrary evaluation of the solution at any point in both time and space.
In what follows, while being misnomer, space and time continuity of the solution designate the continuity of the spatial and temporal domain of definition of the solution, and the continuity of the solution itself.
2

---

Published as a conference paper at ICLR 2024
Continuous solutions for PDEs date back to the early days of deep learning Dissanayake & Phan- Shin-Thien/se 1994 Lagaris et al. 1998 Psichogios & Ungar 1992) and have recently experienced neural networks represent the solution directly as a neural network and train the model to minimize a residual loss derived from the PDE. They are mesh-free, which alleviates the need for complex adaptive mesh refinement tech- niques (mandatory in finite volume methods), and have been successfully applied to broad range of
resurgence of interest Raissi al. (2019]2017). Physics-infonee
physical problems (Lu et al., 2021; Misyris al., 2020 Zoboli et al.) 2022 Kissas et al. 2020; Yang et al., 2019 Cai et al. 2021 with growing community proposing architecture designs specifically
tailored for PDEs (Sitzmann et al. 2020 Fathony et al. (02T) as well as new training methods
(Zeng et al. 2023 Finzi et al. [023] de Avila Belbute-Peree & Kolter 2023
Yet, these models
are also known to be difficult to train efficiently (Krishnapriyan et al.
2021 Wang et al.) 2022).
Recently, neural operators have attempted to learn a mapping between function space, leveraging
kernels in Fourier space (Li et al. 2020b) (FNO) or graphs (Li et al., 2020a) (GNO) to learn the correspondence from the mntiall condition to the solution at fixed horizon. While some operator learning frameworks can theoretically generalize to unseen initial conditions and arbitrary locations, must consider the practical limitations of existing baselines. For instance, FNO requires a static cartesian grid and cannot be directly evaluated outside the training grid. Similarly, GNO can handle arbitrary meshes in theory, but still has limitations in evaluating points outside the training grid and Li et al. (02I)) variant can only be queried at fixed time increments. DeepONet (Lu et al. 2019) can handle free sampling time and space but also constrained to a static observation grid.
Continuous and generalizable solvers- represent a significant challenge. Few models satisfy all these conditions. MP-PDE (Brandstetter et al. 2022) can handle free-form grids but cannot general- ize to different resolutions between train and test, and performs auto-regressive temporal forecasting. Closer to our work, MAgNet (Boussif et al.] 2022) proposes to interpolate the observation graph in latent space to new query points before forecasting the solution using graph neural networks. How- ever, they assume prior knowledge of the evaluation mesh and the new query points, use nearest neighbor interpolation instead of trained attention and struggle to generalize to finer grids during test time. In Hua et al. (2022), the auto-regressive McshGraphNet (Pfaff et al. (0200) is combined with Orthogonal Spline Collocation to allow for arbitrary spatial queries. Finally, DINo (Yin al.) 2022) proposes a mesh-free, space-time continuous model to address PDE solving. The model uses context adaptation techniques to dynamically adapt the output of implicit neural representation forward in time. DINo assumes the existence of a latent ODE modeling the temporal evolution of the context vector and learns it Neural ODE, In contrast, our method differs from DINo as our model is based on physics forecasting in an auto-regressive manner. We achieve space and time continuity through learned dynamical attention transformer capable of handling arbitrary locations and points in time. Our design choices allow for generalization on new spatial and temporal loca- tions, ie. not limited to discrete time steps, and new initial conditions while being trainable from sparse observations
3 CONTINUOUS SOLUTIONS FROM SPARSE OBSERVATIONS
Consider dynamical system following Partial Differential Equation (PDE) defined for all (æ, t) e Q x [0, T7, with T positive constant:
s(xx t) = =f (s,x, t) V(z, t) E f [0,21,
(1)
s(x,00 0) = 8o(x) Vo E Ω s(x,t) = t( V(a,t) E 0n x [0,7]
where the state lies in an invariant set s E S.f S an unknown operator, so + R is the initial condition (IC) and 8 an x [0, + RT the boundary condition. In what follows, we consider trajectories with shared boundary conditions, hence we omit from the notation for readability. In practice, the operator f is unknown, and we assume access to set D of K discrete trajectories from different ICs, st sampled at sparse and scattered locations in time and space. Formally, we introduce two finite sets X' C Ω of fixed positions and fixed regularly sampled times T at sampling rate A Let S(5p, xr 1) be the solution of this PDE from IC 8g, the dataset given
as: D == { (s X,T) k [1, [1. Our task is formulated as:
Given D, new initial condition So E S, and a query (æ, t) E S x equation Mat the queried location and from the given IC, that is S(su, x, t).
10, T1, find the solution of
"Code will he made public. Project page: https: /ttpp:///ootti nuous -pde. iireegittuub.11 /
3

---

Published as a conference paper at ICLR 2024
Latent S Space
Piesocal Staroe
((x,4)
5.f coM
a(a,5)
Prysarcal Space
Taseent Space
Stytem
Time
Plysical Space (Outpaa prodfcctios
(0)
2
Sd [n]
Figure I: Model overview We achieve space and time continuous simulations of physics systems by formulating the task as a double observation problem. System 1 is discrete dynamical model used to compute a sequence of latent anchor states zd auto-regressively, and System is used to design state estimator dig retrieving the dense physical state at arbitrary locations (æ, ().
Note that this task involves generalization to new ICs, as well as estimation to unseen spatial loca- tions within Ω and unseen time instants within [0, 7T. We do not explicitly require extrapolation to instants T, although comes as side benefit of our approach up to some extent.
3.1 THE DOUBLE OBSERVATION PROBLEM
The task implies extracting regularities from weakly informative physical variables that are sparsely measured in space and time, since X' and T contain very few elements. Consequently, the possi- bility to forecast their trajectories from off-the-shelf auto-regressvve methods is very unlikely (as confirmed experimentally). To tackle this challenge, we propose an approach accounting for the fact that the phenomenon is not directly observable from the sparse trajectories, but can be deduced from a richer latent state-space in which the dynamics is markovian. We introduce two linked dynami- cal models lifting sparse observations to dense trajectories guided by observability considerations, namely
System 1: za[n+]] fi (zal(n [m]) s2
= h1 (zaan) System 2: {
s(x,t) = =f, V(x, t)e0 xxx,t) = h2 (s,
(2) where for all N. we note sain] s(X, nd) the sparse observation at some instant n (the sampling rate A is not necessarily equal to the sampling rate A used for data acquisition, which we will exploit during training to improve generalization. This will be detailed later).
System 1 1- is a discrete-time dynamical system where the available measurements s2(n] are con- sidered as partial observations of a latent state variable za[n].. We aim to derive output predictor from System I to forecast trajectories of sparse observations auto-regressively from the sparse IC. As mentioned earlier, sparse observations are unlikely to be sufficient to perform predictions, hence we introduce a richer latent state variable zd in which the dynamics is truly markovian, and obser- vations sa[r] are seen as measurements of the state zd using the function fij.
System 2 - is a continuous-time dynamical system describing the evolution of the to-be-predicted dense trajectory S(Sp, xr t). introduces continuous observations 2(x f) such that =xX,,, n.)) za[n]. The insight is that the state representation =a[n] obtained from System is designed to contain sufficient information to predict sa[n], but not necessarily to predict the dense state. Formally, zd represents solely the observable part of the state, in the sense of control theory.
At inference time, we forecast at query location (ar, t) with a 2-step algorithm: (Step-1) System 1 is used as output predictor from the sparse IC (0). and computes a sequence 2(0), = ..... which we refer to as "anchor states". This sequence allows the dynamics to Markovian, provides sufficient information for the second state estimation step and holds information to predict the sparse observations, allowing supervision during training. (Step-2) We derive a state observer from System 2 leveraging the anchor states arbitrary location in space and time (see figure I]. Importantly, for a given IC, the anchor states are computed only once and reused within System 2 to estimate the solution at different points.
the whole time domain to estimate the dense solution at an
4

---

Published as a conference paper at ICLR 2024
3.2 THEORETICAL ANALYSIS
In this section, we introduce theoretical results supporting the use of Systems and 2. In particular, show that using System 1 to forecast the sparse observations in latent space za rather than directly operating in the physical space leads to smaller upper bounds on the prediction error. Then, we show the existence of a state estimator from System 2 and compute an upper bound on the estimation error depending on the length of the sequence of anchor states.
Step 1- consists of computing the sequence of anchor states guided by an output prediction task of the sparse observations. As classically done, we introduce an encoder (formally, a state observer) e( ((sa(0)) z [0] coupled System to project the sparse IC into a latent space za. Following System 1, we compute the anchor states Zat auto-regressively (with fi) in the latent space. The sparse observations are extracted from 2 using f1 In comparison, existing baselines (Pfaft et al.) 2020 -oohallz et al., 2020 Stachenfeld al. 202I) maintain the state in the physical space and discard the intermediate latent representation between iterations. Formally, let us consider approximations fú, fi, è (in practice realized as deep networks trained from data D) of fi, hı and e and compare the prediction algorithm for the classic auto-regressive (AR) approach and ours
Classic AR: 8 [x] == (hi o fi ē)" (s)(0)
Ours:
Sa[n] == hi fi " o é (sx(00)
(3)
Classical AR approaches re-project the latent state into the physical space at each step and repeat Our method encodes the sparse IC, advances the system in the latent space, and decodes toward the physical space at the end. A similar approach has also been explored in Wu et al. (2022); Kochkov et al. (2020). albeit in different contexts, without theoretical analysis.
sencode-process-decee
Proposition 1 Consider a dynamical system of the form of System I and assume the existence of state observer along with approximations f1, hi, with Lipschitz constants L,,, Ln and Le respec- tively such that Ln L Le + 1. If there exist 8y, 5n, 8e E R such that V(z 8) E R== x RR,
[z(()) - (z)) < og,
le(s) - ε(s) < 8e |, then for all integer n > O, with Sa[n] and 87 [m] as in equation
|/1 (z)) & dn.
(4)
for the Euclidean norm
85 [.
[8a[r]] 8a|mi]1 < 8n + Ln
+ 1.750,)
(5)
T.g -
6 [x]] < -1 --1 with 8 = Ön + Lno, + Lyd, and L = Ln Lygl
sa|n]] -
(6)
Proof: See appendix B
This result shows that that falling back to the physical space at each time step degrades the upper bound of the prediction error. Indeed, if T. < 1, the upper bound converges trivially to zero when n increases, and hence can be ignored. Otherwise, the upper bound for the classic AR scheme appears to be more sensitive to approximation errors o,, 6y and & compared to our approach (for a formal comparison, see appendix C]. Intuitively it means that information is lost in the observation space, which thus needs to be re- re-estimated at each iteration when using the classic AR scheme. By maintaining a state variable the latent space, we allow this information flow readily between each step of the simulator (see blue frame in figure D..
Step The state estimator builds upon System 2 and relies on the set of anchor states from the previous step to estimate the dense physical state at arbitrary locations in space and time. Formally, we look for a function 4/g leveraging the sequence of anchor states zu [0]. -.)..... =x[g] (simulated from the sparse IC [a) to retrieve the dense solution In what follows, we show that (1) such a function
vi, exists and (2) we compute upper bound on the estimation error depending on the length of the sequence. To do so, consider the functional which outputs the anchor states from any IC 80 E
(p()))= hug (x)) by (S(s,, X, A)) ... (S X, X))A)) ]-[ 20
zapp]
(7)
In practice the ground truths =n[n]] are not perfectly known, as they are obtained from data-driven output predictor (step 1) using the sparse IC. Inspired from Janny al. (0022)), we state:
Since the simulation is conducted up to T, and considering the time step A, in practice q & 1
5

---

Published as a conference paper at ICLR 2024
Proposition 2 Consider dynamical system defined by System 2 and equation Assume that
AL. fo is Lipschitz with constant Le.
A2. there exists p > O and a strictly increasing function a such that Vsa. 86 S and Vq b p [0,(sa) 0y((6)] > a(g)]sa sols
(8)
where| Is is an appropriate norm for S.
Then, 4q 2 p, there exists Vig such that, for (xx t) E 12x [0, I'] and on such that Za(n] =a[n] for all n & q.
Vig (z [0), ... zu(g], xz Ss0, 1. S(00 1) 1h , 24, 1)/ < <2x(1) (oo1))"(ooogle"-
.M zalu),
(10)
where 601==[60 8g].
Proof: See appendix D Assumption A2. states that the longer observe two trajectories from different ICs, the casier it will be to distinguish them, ruling out systems collapsing to the same state. Such systems are uncommon since forecasting their trajectory becomes trivial after some time. This assumption related to finite-horizon observability in control theory, a property of dynamical systems guaranteeing that the (markovian) state can be retrieved given a finite number p of past observations. Equation g is associated with injectivity of Oq, hence the existence of a left inverse mapping the sequence of anchor states to the IC su-
Proposition 2 highlights trade-off on the performance of eiy. On one hand. longer sequences of anchor states are harder to predict, leading negatively. On the other hand, longer sequences hold more information that can still be leveraged by Vg to improve its estimation, represented by aa)) -1 in equation 10 contrast to competing baselines or conventional interpolation algorithms, our approach takes this trade-off into account, by explicitly leveraging the sequence to estimate the dense solution, as will be discussed below.
larger 1ōo which impacts the state estimator vi
Discussion and related work the competing baselines can be analyzed using our setup, yet in weaker configuration. For instance, one can see Step 2 as an interpolation process, and replace with a conventional interpolation algorithm, which typically relies spatial neighbors only. Our method not only exploits spatial neighborhoods but also leverages temporal data, improving the performance, shown in proposition and empirically corroborated in Section 4
MAgNet (Boussif et al., 2022) uses a reversed citrpolate-forecaa scheme compared to ours. The IC sa is interpolated right from the start to estimate su (corresponding to Step 2, with q=1), and then simulated with an auto-regressive model in the physical space (with the classic AR scheme). Propositions 1 and [ show that the upper bounds the estimation and prediction error are higher than ours. Moreover, if the number of query points exceeds the number known points ([]}||)) the input of the regressveee solver is filled with noisy interpolations, which impacts performance.
DINo (Yin et al. 2022) is a very different approach leveraging spatial implicit neural representation modulated by a context vector, whose dynamics is modeled via a learned ODE. This approach is radically different than ours and arguably involves stronger hypotheses, such as the existence of a learnable ODE modeling the dynamics of suitable weight modulation vector. method relies on arguably more sound assumptions, i.e. the existence of an observable discrete dynamics explaining the sparse observation, and the finite-time observability of System 2.
contrast, our
3.3 IMPLEMENTATION
The implementation follows the algorithm described in the previous section: (Step-1) rolls out pre- dictions of anchor states from the IC, (Step-2) estimates the state at the query position from these anchor states. The encoder é from Step 1 is a multi-layer perceptron (MLP) which takes as in- put the sparse IC sa(0] and the positions X' and outputs latent state variable =,[0] structured as a graph, with edges computed with a Delaunay triangulation. Hence, each anchor is a graph =a[n] = { but we will omit index i over graph nodes in what follows if not required for understanding.
We model fi as a multi-layer Graph Neural Network (GNN) (Battaglia et al. 2016). The anchor states [n[n] are defined at fixed time steps nA,, which might not match A used in the data T. We

---

Published as a conference paper at ICLR 2024
found it beneficial to choose A=kxA" with k>1eN such that the model can be queried during training on time points E T that do not match exactly with every time-steps in ==[0], 2 [1] but rather subset of them, hence encouraging generalization to unseen time. The observation function hi is MLP applied on the vector at node level the graph z.t
The state estimator vi, is decomposed into Transformer model (Vaswani et al. 2017) coupled to a recurrent neural network to provide an estimate at query spatio-temporal query position (xc, (). First, through cross-attes -tossntttnntoon we translate the set of anchor states =a[nz] (one embedding per graph node i and per instant n) into n., which we denote 2nd (æx, t) (one embedding per instant n)). Following advances in geometric mappings in computer vision (Saha et al.) (022). we use multi-head cross-attention to query from coordinates (ær, t) to Keys corresponding to the nodes in each graph anchor state za(n), Vn:
of estimates of the continuous variable 2 2(x t) conditioned at the instant
Znd (x. t) = fmna (Q- t), (=V=<zamm) + Co, nA)). // attention over nodes i (11)
where Q, K, V are, respectively, Query, Key and Value inputs to the cross-attention layer fiha Vaswani et al. 2017) and Sw Fourier positional encoding with a learned frequency parameter w. Finally, we leverage a state observer to estimate the dense solution at the query point from the time. This is achieved with a Gated Recurrent Unit (GRU) Cho et al. (2114) maintaining a hidden state us(n).
sequence of conditioned anchor variables,
u/n] ggu (u (uIn-1]. znA (x. t)), S(s). x. D (u(g)),
(12)
which shares similarities with conventional state-observer designs in control theory (Bernard et al., 2022). Finally, an MLP D maps the final GRU hidden state to the desired output, that is, the value of the solution at the desired spatio-temporal coordinate (xx, t). See appendix E for details.
3.4 TRAINING
Generalization to new input locations during training is promoted by creating artificial generalization situations using sub-sampling techniques of the sparse sets X' and T.
Artificial generalization The anchor states za[ri] are computed at time rate larger than the available rate A This creates situations during training where the state estimator 3y does not have access to a latent state perfectly matching with the queried time. We propose a similar trick promote spatial generalization. At each iteration, we sub-sample the (already sparse) IC sa[0] randomly to obtain [2 defined on a subset of X. We then compute the anchor states žd using System On the other hand, the query points are selected in the larger set X. Consequently, System 2 is exposed to positions that do not always match with the ones in =a(r]. Note that the complete domain of definition 92 [O, remains unseen during training.
Training objective To reduce training time, we randomly sample M query points (Im. Tm) in X' T at each iteration. with probability proportional the previous error of the model this point since its last selection (see appendix and we minimize the loss
Lomnans
Laraanic
IT/Al
S( Fm, Tm) - 2h (Eal0/a), x )
Sa[n] - hi (Zalmi)
L
+
(13)
with _n[n]] f x à (5)(0)). Ceontimoous supervises the model end-to-end, and Loynamics trains the latent anchor states za to predict the sparse observations from the IC.
4 EXPERIMENTAL RESULTS
Experimental setup X'x results from sub-sampling 2 7] with different rates to control the difficulty of the task. We evaluate on three highly challenging datasets (details in appendix (F) Navier (Yin et al.) 2022 Stokes 2009) simulates the vorticity of viscous, incompressible flow driven by a sinusoidal force acting on a square domain with periodic boundary conditions. Shallow Water (Yin et al. 2022 Galewsky et al. 2004) studies the velocity of shallow waters evolving the tangent surface of a 3D sphere. Eagle (Janny et al. 2023) is a challenging dataset of turbulent airflow generated by moving drone in 2D environment with many different scene geometries.
We evaluate our model against three baselines representing the state-of the-trr in continuous simu- lations. Interpolated MeshGraphNet (MGN) (Praff al. 2020) is standard multi-layered GNN
7

---

Published as a conference paper at ICLR 2024
Navier
Shallow Water
Eagle
High
Mid
Low
High
Mid
Low
High
Low
DINo
In-X
I,557
1.130
1.878
0.1750
0.1814
0.2733
287.3
302.7
(Yin ct al. 2022)
Ext-X
1.600
1.253
5.493
4.638
13.40
21.55
381.7
489.6
Interp. MGN
In-X
1.913
0.9969
0.6012
0.3663
0.2835
0.7309
64.44
83,58
(Pfaff et al. 2020)
Ext-X'
2.694
4.784
14.80
1.744
4.221
8.187
173.4
241.5
In-X
n/a
n/a
n/a
Time Oracle (n.c)
Ext-X
0.851
4.204
15.63
1.617
4.327
8.522
147.0
221.2
MAgNet
In-X
18.17
6.047
8.679
0.3196
0.3358
0.4292
99.79
124.5
(Boussif et al. 2022)
Ext-X
35.73
26.24
57.21
10.21
23.20
30.55
194.3
260.7
In-X
0.1989
0.2136
0.2446
0.2940
0.3139
0.2700
70.02
78,83
Ours
Ext-X'
0.2029
0.2463
0.5601
0.4493
1.051
2.800
90.88
117.2
Table 1: Space Continuity we evaluate the spatial interpolation power of our method vs. the baselines and standard interpolation techniques. We vary the number of available measurement points in the data for training from High (25% of simulation grid), Middle (10%), and Low (5%) amount of points and show that our model outperforms the baselines. Evaluation is conducted over 20 frames in the future (10 for Eagle) and we report the MSE to the ground truth solution (x10-9).
used reto-seegrssivvll and extended spatiotemporal continuity using physics-agnostie interpola- tion. MAgNet (Boussif et al. 2022) interpolates the IC at the query position in latent space before using MGN. The original implementation assumes knowledge of the target graph during training, including new queries. When used for superresolution, the authors kept the ratio between the amount of new query points and available points constant. Hence, while MAgNet is queried at unseen loca- tions, it also benefits from more information. In our setup, the model is exposed to fixed number of points but does not receive more samples during evaluation. This makes our problem more challeng- ing than the one addressed in Boussif et al. (2022). DINo (Yin et al.) 2022) models the solution as an Implicit Neural Representation (INR) s(x at) where the spatial coordinates x are fed to a MFN (Fathony et al. 202T) and a; is context vector modulating the weights of the INR. The dynamics of a is modeled with a Neural-ODE, where the dynamics is an MLP. We share common objectives with DINo and take inspiration from their evaluation tasks yet in more challenging setup. Details of the baselines are in appendix F We highlight a caveat on MAgNet: the model can handle a lim- ited amount of new queries, roughly equal the number of observed points. Our task requires the solution at up to 20 times more queries than available points. In this situation, the graph in MaGNet dominated by noisy states from interpolation, and the auto-regressive forecaster performs poorly. During evaluation, we found it beneficial to split the queries into chunks of 10 nodes and to apply the model several times. This strongly improves the performance at the cost of increased runtime.
Space Continuity - Table compares the spatial interpolation power of our method versus several baselines. The MSE values computed on the training domain (In-X'=X') and outside (Ext-X'=(\) ( show that our method offers the best performance, especially for the Ext-domain task, which is our aim. To ablate dynamics and evaluate the impact of trained interpolations, we also report the pre- dictions of Time Oracle which uses sparse ground truth values at all time steps and interpolates (bicubic) spatially. This allows us to assess whether the method is doing better than a simple ax- iomatic interpolation. While MGN offers competitive in-domain predictions, the cubic interpolation fails to extrapolate reliably on unseen points. This can be seen in the In/Ext gap for Interpolated MGN which very close the Time Oracle error. MaGNet, which builds on similar framework, is hindered by the larger amount of unobserved data in the input mesh. At test time, the same num- ber of initial condition points are provided but the method interpolates substantially more points. DINo achieves a very low In/Ext gap. yet fails on highly (5%) down-sampled tasks. One of the key difference with DINo that the dynamics relies on an internal ODE for the temporal evolution of a modulation vector. In contrast, our model uses an explicit auto-regressive backbone, and time forecasting is handled in an arguably more meaningful space, which we conjecture to be the reason why we achieve better results (see fig. 5 in the appendix).
Time Continuity - is a step forward in difficulty, as the model needs to interpolate not only to unseen spatial locations (datasets are undersampled at 25%) but also on intermediate timesteps (Ext- T. Table 2) All models perform well on Shallow Water, which is relatively easy. Both DINo and MAgNet leverage discrete integration scheme (Euler for MAgNet and RK4 for DINo) allowing querying the model between timesteps seen at training. These schemes struggle to capture the data dependencies effectively and therefore the methods fail Navier (see also Figure 6 for qualitative
8

---

Published as a conference paper at ICLR 2024
to +2
+
to
to +-1
10
13
DCN
20
15
MAg
Figure 2: Results Eagle Per point error of the flow prediction Eagle example in the Low
spatial down-sampling scenario. Our model exhibits lower errors as also shown in Tables and
Navier
Shallow Water
Eagle
1/1
1/2
1/4
1/1
1/2
1/4
1/1
1/2
1/4
DINo
In-T
1.590
36.31
46.02
3.551
6.005
6.249
444.5
447.1
448.6
YYn et al. 2022)
Ext-T
n/a
39.42
54.72
n/a
6.015
6.265
n/a
479.4
470.7
Interp. MGN
In-T
2,506
4.834
12.77
L.488
1.289
1.333
203.4
210.4
263.3
(Pafff et al. 2020)
Ext-T
n/a
5.922
36.43
n/a
1.287
1.355
n/a
209.8
263.8
In-7
n/a
n/a
n/a
Spatial Oracle (n.c)
Ext-T
n/a
1.296
28.58
n/u 0.003
0.119
n/a
29.46
54.53
MAgNet
In-7
31.51
135.0
243.9
7.804
6.433
1.884
227.9
220.3
225.8
(Boussif et al. 2022
Ext-T
n/a
142.8
255.5
n/a
6.291
1.947
229.8
230.6
In-T
0.2019
0.1964
0,4062
0.4115
0.4278
0.4549
108.0
106.1
278.6
Ours
Ext-T
n/a
0.2138
11.36
n/a
0.4326
0.4802
n/a
119.9
306.9
evaluate the time interpolation power of our method vs. the base- lines. Models are trained and evaluated with 25% of Q and with different temporal resolutions (full, half, and quarter of the original). The Spatial Oracle (not comparable!) uses the exact solution at every point in space, and performs temporal interpolation. Evaluation is conducted over 20 frames the future (10 for Eagle) and we report MSE compared the ground truth solution (x 10)
Table 2:
Time Continuity
results). Eagle is particularly challenging, the main source of error being the spatial interpolation, as can be seen in Figure our method yields lower errors in flow estimation.
Many more experiments are available in appendix G We study the impact of key design choices, artificial generalization, and dynamical loss. We show qualitative results on time interpolation, time extrapolation on the Navier dataset. We explore generalization to dif- ferent grids. We provide more empirical evidence of the soundness of Step in an ablation study (including comparison with attentive neural process Kim et al. (0118). attention- based structure somehow close to ours), and observe attention maps on several examples. We show that our state estimator goes beyond local interpolation, polation algorithms would do. Finally, we also measure the computational burden of the discussed methods and show that our approach more efficient.
conventional inter-
5 CONCLUSION
We exploit a double dynamical system formulation for simulating physical phenomena at arbitrary locations in time and space. Our approach comes with theoretical guarantees existence and ac- curacy without knowledge of the underlying PDE. Furthermore, our method generalizes to unseen initial conditions and reaches excellent performances outperforming existing methods. Potential applications our model goes beyond fluid dynamics and can be applied to various PDE-based problem. Yet. our approach relies several hypotheses such regular time sampling and observ- ability. Finally, for known and well-studied phenomena, it would interesting to add physics priors in the system, nontrivial extension that we leave for future work.

---

Published as a conference paper at ICLR 2024
Reproducibility the detailed model architecture is described in the appendix. For the sake of reproducibility, in the case of acceptance, we will provide the source code for training and evaluating model, as well as trained model weights. For training, up the codebase, including installing external dependencies, pre-trained models, and pre-selected hyperparameter configuration. For the evaluation, the code will include evaluation metries directly comparable to the paper's s results.
will provide instructions for setting
Ethies statement While our simulation tool is unlikely yield unethical results, we are mindful of potential negative applications of improving fluid dynamics simulations, particularly in military contexts. Additionally, we strive to minimizing the carbon footprint associated with our training processes.
6 ACKNOWLEDGEMENTS
We recognize support through French grants "Delicio" (ANR-19-CE23-0000 of call CE23 "Intelli- gence Artificielle" and "Remember" (ANR-20-CHIA0015)) of call "Chaires IA hors centres". This work was performed using HPC resources from GENCI- IDRIS (Grant 2023-AD010614014)
REFERENCES
Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al. Interaction networks for learning about objects, relations and physics. Neural Information Processing Systems, 2016.
Pauline Bernard, Vincent Andricu, and Daniele Astolfi. Observer design for continuous dy- namical systems. Annual Reviews in Control, 2022.
Emmanuel De Bézenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical processes: Incorporating prior scientific knowledge. Journal of Statistical Mechanics: Theory and Experi- ment, 2019.
Oussama Boussif, Yoshua Bengio, Loubna Benabbou, and Dan Assouline. Magnet: Mesh agnostic neural pde solver. In Neural Information Processing Systems, 2022.
Johannes Brandstetter, Daniel E. Worrall, and Max Welling, Message passing neural PDE solvers. In International Conference on Learning Representations, 2022.
Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physics- informed neural networks (pinns) for fluid mechanies: A review. Acta Mechanica Sinica, 2021.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Neural Information Processing Systems, 2018.
Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Hol- ger Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoder-decoder for statistical machine translation. arXiv preprint, 2014.
Filipe de Avila Belbute-Peres and J Zico Kolter. Simple initialization and parametrization of sinu- soidal networks via their kernel bandwidth. In International Conference on Learning Represen- tations. 2023.
MWMG Dissanayake and Nhan Phan-Thien. Neural-network-based approximations for solving partial differential equations. Communications in Numerical Methods in Engineering, 1994.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. Neural Information Processing Systems, 2019.
Rizal Fathony, Anit Kumar Sahu, Devin Willmott, and J Zico Kolter. Multiplicative filter networks. In International Conference on Learning Representations, 2021.
Mare Anton Finzi, Andres Potapezynski, Matthew Choptuik, and Andrew Gordon Wilson. A sta- ble and scalable method for solving initial value pdes with neural networks. In International Conference Learning Representations, 2023.
10

---

Published as a conference paper at ICLR 2024
Joseph Galewsky, Richard K. Scott, and Lorenzo M. Polvani. An initial-value problem for testing numerical models of the global shallow-war watter equations. Tellus A: Dynamic Meteorology and Oceanography, 2004.
Vincent Le Guen and Nicolas Thome. Disentangling physical dynamics from unknown factors for unsupervised video prediction. In Conference on Computer Vision and Pattern Recognition, 2020.
Xu Han, Han Gao, Tobias Pfaff, Jian-Xun Wang, and Liping Liu. Predicting physics in mesh- reduced space with temporal attention. In International Conference on Learning Representations, 2021.
Chuanbo Hua, Federico Berto, Michael Poli, Stefano Massaroli, and Jinkyoo Park. Efficient con- tinuous spatio-temporal simulation with graph spline networks. In Internation Conference on Machine Learning (AI for Science Workshop), 2022.
Steeven Janny, Fabien Baradel, Natalia Neverova, Madiha Nadri, Greg Mori, and Christian Wolf. Filtered-cophy: Unsupervised learning of counterfactual physics in pixel space. In International Conference on Learning Representation, 2022a.
Steeven Janny, Quentin Possamai, Laurent Bako, Christian Wolf, and Madiha Nadri. Learning reduced nonlinear state-space models: an output-error based canonical approach. In Conference on Decision and Control, 2022b.
Steeven Janny, Aurélien Beneteau, Nicolas Thome, Madiha Nadri, Julie Digne, and Christian Wolf. Eagle: Large-scale learning of turbulent fluid dynamics with mesh transformers. In International Conference on Learning Representation, 2023.
Hyunjik Kim, Andriy Mnih, Jonathan Schwarz, Marta Gamelo, Ali Eslami, Dan Rosenbaum, Oriol Vinyals, and Yee Whye Teh. Attentive neural processes. In International Conference on Learning Representations, 2018.
Georgios Kissas, Yibo Yang, Eileen Hwuang, Walter Witschey, John A. Detre, and Paris Perdikaris. Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive flow mri data using in Applied Mechanics and Engineering, 2020.
neural networks. Computer Methods
Dmitrii Kochkov, Alvaro Sanchez-Gonzalez, and Peter Battaglia. Learning latent field dynamics of pdes. In Third Workshop on Machine Learning and the Physical Sciences (NeurIPS 2020), 2020.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Char- acterizing possible failure modes physics-informed neural networks. Neural Information Pro- cessing Systems, 2021.
Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. Transactions Neural Networks, 1998.
Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba. Learning par- ticle dynamics for manipulating rigid bodies, deformable objects, and fluids. In International Conference on Learning Representations, 2018.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhat- tacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differ- ential equations. In Neural Information Processing Systems, 2020a.
Zongyi Li. Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al. Fourier neural operator for parametric partial differential equa- tions. In International Conference on Learning Representations, 2020b.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu. Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Markov neural operators for learning chaotic systems. arXiv preprint, 2021.
11

---

Published as a conference paper at ICLR 2024
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for iden- tifying differential equations based on the universal approximation theorem of operators. arXiv preprint, 2019.
Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven G Johnson. Physics-informed neural networks with hard constraints for inverse design. Journal on Scientific Computing. 2021.
George Misyris, Andreas Venzke, and Spyros Chatzivasileiadis. Physics-informed neural networks for power systems. In Power & Energy Society General Meeting, 2020.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzaler, and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2020.
Dimitris C Psichogios and Lyle H Ungar. A hybrid neural aetwork-first principles approach to process modeling. AIChE Journal, 1992.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physies informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint, 2017.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 2019.
Prajit Ramachandran, Barret Zoph, and Quoc V Le. Searching for activation functions. arXiv preprint, 2017.
Avishkar Saha, Oscar Mendez, Chris Russell, and Richard Bowden. Translating images into maps. In International Conference on Robotics and Automation, 2022.
Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Confer- ence on Machine Learning. 2020.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Neural Information Processing Systems, 2020.
Kim Stachenfeld, Drummond Buschman Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Sonzhez- Learned simulators for turbulence. In International Conference on Learning Representations, 2021.
George Gabriel Stokes. On the Effect of the Internal Friction of Fluids on the Motion of Pendulums. Cambridge University Press, 2009.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Neural Information Processing Systems, 2017.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail train: A neural tangent kernel perspective. Journal of Computational Physics, 2022.
Tailin Wu, Takashi Maruyama, and Jure Leskovee. Learning to accelerate partial differential equa- tions via latent global evolution. Advances in Neural Information Processing Systems, 2022.
X. I. A. Yang, S. Zafar, J.-X. Wang, and H. Xiao. Predictive large-eddy-simulatin wall modeling via physics-infee inforree neural networks. Physical Review Fluids, 2019.
Yuan Yin, Matthieu Kirchmeyer, ovan-Yves Franceschi, Alain Rakotomamonjy, et al. Continuous pde dynamics forecasting with implicit neural representations. In International Conference on Learning Representations, 2022.
Qi Zeng, Yash Kothari, Spencer H. Bryngelson, and Florian Schäfer. Competitive physics informed networks. In International Conference on Learning Representations, 2023.
Samuele Zoboli, Steeven Janny, and Mattia Giaccagli. Deep learning- based output tracking via regulation and contraction theory. In International Federation of Automatic Control, 2022.
12

---

Published as a conference paper at ICLR 2024
A WEBSITE AND INTERACTIVE ONLINE VISUALIZATION
An anonymous website has been created where results can be visualized with an online interactive tool, which allows to choose time steps interactively with the mouse, and in the case of the Shallow Water dataset, also the orientation of the spherical data:
https:://cottiiuuuuspddd.e github. io/
B PROOF OF PROPOSITION
The proof proceeds by successive majorations and triangular inequalities, For sake of clarity, and only in this proof we omit the d subscript and write s[n] and x/n]] for sa[n] and =a[n],, respectively.
We start with 8(m] == fú A m é ((s())) Thus for any integer n we have
s|n]] - 8[r]]] //1 (z[n²) hi (2|ml)]-
(14)
Using Lipschitz property and then
[s|n] - [[m]] < |/i (=[ni) - hi (z[n]) + 1 (=[ml) hi (2[n))
(15)
8 L [n] _|r]]|
Noticing that one can rewrite 2 [/] as [[n] = f o {s(0)) Since z[n] = f1 (z[0)) and using a similar decomposition as for 15) one gets:
|||n| ||||||
< 8, Σ
LÝ + -3]= =|0]|-
(16)
k=U
Hence, from equation 15 and using 2(0] = e(s(0)) and 2|0] = é (s(())) we have
Li -
[ s[n]] & 5n + Ta 6 1,, 1 + 1.75e
(17)
We now move to the classic auto-regressive case, i.e. sar ["[[n] (ii fa 2) (s]0))
in]
s" [n]]
< h1 (=[nl) hi (z(nl) + [2([l) hn (1))
(18)
< 5 + In =|m] ==[m]]
6y + Ls cs/m - 11] 6(8" -111
< 8, + Ln
< 8, + Ln
6, + Ls (5e + Le s//n - 1] - 8 [n 111)
n-2
<
Σ
L' + |s|[1] - s.H [1]]
1=0
with 5 6n + Lnd5 + Ln 5, and I, I.,, I. Lec. Moreover,
s{1]] - [1]
= |// (z (1] - f (1))
(19)
< 6n La [ - [1]
(20)
< 6n + La + Ly|x[0] -- (0])
(21)
< 5
(22)
Putting it all together, we get equation 6
1
[s[n]]
8 [n]]
<
(23)
T,
Finally, equation 17 and equation 23 conclude the proof.
13

---

Published as a conference paper at ICLR 2024
C COMPARISON OF UPPER BOUNDS IN PROPOSITION I
We start by formulating equation 5 and equation
6 under a comparable form
L
1 + La 5 1) Lg
[ 8a[na] < 8 + Ln Ln
= Ki
(24)
1 1ー-1 L L
[sa[n]
22 [r]]
8 K2
(25)
Now we consider two cases depending on the Lipschitz constants of the problem. namely Lns Lgs and Le. Finst, consider the case where the Lipschitz constants are very large (i.e. Lis Lig, Le >> 1). In that case, the upper bounds can be approached by
K1 Ln + LnLjde
(26)
K2 se ón d
Lnd Ln LE + -
(27)
Hence, K₂ > K1 (we highlighted the difference between both terms in the previous equation. Now consider the case where the Lipschitz constants are very small (i.e. Lxs Ly. Le < 1). Recall that this case corresponds to a trivial prediction task since any trajectory of System will converge to unique state. Again, the upper bounds can be approached by
K1 ~
(28)
K₂ L5
(29)
In this trivial case, the upper bound on the prediction error using our method is a combination of the approximation errors from each function. On the other hand, using the classic AR scheme implies a larger error, since the model accumulates approximations at each time step not only from the dynamics but also from the observation function and the encoder.
D PROOF OF PROPOSITION 2
The proof follows the lines of Janny et al. 2022b). The existence of bh is granted by the observ-
ability assumption. Indeed assumption A2. states that for all p. Og is injective in S. Hence, such that Vs' E S
exists a inverse mapping O Og (O, =s'
(30)
Let za ==[0|g]=[<ad]] zalgi]. Hence, one can build bhy using the dynamics of the system for all Ω: Yso E S, S(50, S(; x,
= vi (za [0l] xx.)
(31)
Now, because of the noise, the disturbed observation 2.[0/q] za[0|g] + 6019 may not belong to (S), where the inverse mapping O is well defined. We solve this by finding the closest "possible" observation.
an = arg min
e'ES [a(0 0,(s')
(32)
s(x, t) S(so x, == v. (2,1000)) x).
(33)
Hence, we have for all s' S
La 00g] - 0g
< |2[00g] -0((s')
(34)
In particular, for s' so and since (,())) za(0 q).
1/ - O, <
1d[0 O,(50)
(35)
In the other hand, from assumption A2. equation
alp)( - sols < [0,(8)) (,
(36)
(0,
2x(0|g]| + [2][ O,(s))
2 200gal
14

---

Published as a conference paper at ICLR 2024
Arbitrary query xXx)
Piysa Spart
Pomnee Pin East
MILA
MHA
Sabriee x queay 4X.
Pre-itmotop
X
draazis teteppaazio
Plyics Space Spanse intial conditiion
Latent Spane
* ahazed weights
Ext-X
Auto-rogrosiive forerasting In latent space
2a[1]
z. [2]
loot
Phyics Space
Physies grounding
Txt.T
Supervised with obecrvations
from
Figure 3: Model overview The model leverages dynamical system (System 1) perform auto- regressive predictions of the dynamics in mesh-structured latent space from sparse imitial condi- tions. It is combined with data-driven state estimator derived from another continuous-ting time dy- namical system (System 2), implemented with multi-head cross-atiention. The attention mechanism queries the intermediate anchor states from the auto-regressive predictor and uses Fourier positional encoding encode the query points (x, r). An additional GRU refines the dynamics after interpo- lation.
Moreover, since f2 Lipschitz
(S0. x) - S(S0, t)ls = =152(S(80 xz ()) - -2((S(0n x. l) Ot
(37)
< LL]S((su x) - S(SB: ac, t)ls
and using the Grönwall inequality
[S(8.. L) S( 1)s & ettallan 80lss
(38)
Finally, combining equation
36 and equation 38
[S( t) S(So x, t)ls & 2a(g) (0000)"'0000(("""
which concludes the proof.
E MODEL DESCRIPTION
In this section, we describe the architecture of our implementation in more detail. The model is illustrated in figure
Step The output predictor derived from System implemented as multi-layer graph neural network inspired from Pfaff et al, (2020); Gonzalzz-Gonaaezz et al. (2020) but without following the standard "encode "encode-process-derr" process-decode"" setup. Let X' {Zo, .... Kx be the set of sub-sampled positions extracted from the known locations (cf. Artificial generalization from section 3.4). The input of the module is the initial condition at the sampled points and the corresponding positions (a,, [2[0 and encoded into a graph-structun
latent space [0 = (z,[0),, e(0)ey)s.s where ==[0]: is a latent node embedding for position æ; and e(0)sş an edge embedding for edge pairs (i, extracted from Delaunay triangulation. The encoder maps the sparse IC node and edge embeddings using two MIPss, Jodge and Snode:
za ==(0]; fmote (s:))) a0;), e e(0)ij Sedge (ar: D.j las 2z51)s
(39)
15

---

Published as a conference paper at ICLR 2024
fuodes and fedges two ReLU-activated MLPs, each consisting of 2 layers with 128 neurons. The initial node and edge features =x[0]; and e(0]s) are represented as 128-dimensional vectors.
The dynamics f is modeled as a multi-layered graph neural network inspired from Pfaff et al. 2020
Sanchez- Gonzalez Sanchez-Gonzalez al. (2020), we therefore add a layer superscript ' to the notation:
s+1
Eig + Bodgo
(=1, zj. eig):
*a(n 1] = fi (za(nl) (xf,et))., such that
eft)
xf + ghane
2 E, Etj
(40)
= e(nliss
= zd(n]is
The GNNs employ two MLPs Foode and Pedges with same dimensions as fodges and fnodes. We compute the sequence of anchor states zd [0]. ... za(g] q) in the latent space by applying fi auto-regressively.
The observation function hi extracts the sparse observations 8a[r]] from the latent state za(nz] and consists of a two-layered MLP with 128 neurons, with Swish activation functions (Ramachandran al. 2017) applied on the node features, i.e. (x:))( Rs lig (z(m]:).
Step - The spatial and temporal domains Ω [00, are normalized, since it tends improve generalization on unseen locations. The state estimator bi takes as input the sequence of latent graph representation =z[0].-.. za|q] and spatiotemporal query sampled in x [0. T]. This query is embedded in a Fourier space using the function Ew which depends on frequency parameter E Idim 1 (initialized uniformly [0. 1)). By concatenating harmonies of this frequency up to some rank, we obtain a resulting embedding of 128 dimensions (if (ar, () exceeds the number of dimensions, cropping is performed to match the target shape).
Ca t) [ cow(wo/ss))))) |n. x), sin skow (nn x), cos(kw,,, +1t), sin +1t), 10, (.....K) ..).
(41)
The continuous variables AnA (ae, 1) conditioned by the anchor states are computed with a multi-head attention Vaswani et al. (001)
2nA (2, t) = funu (Q= (x, t). K-V
(42)
where fnha is defined as
qi A(Q,K,V)
Q2 Q Q+q1,
(43)
Q3 B(gz);
out q3 +
Here, -(-,..) refers to the multi-head attention mechanism described in Vaswani et al. with four attention heads, and B()) represents a single-layer multi-layer perceptron activated by the rectified linear unit (ReLU) function. We do not use layer normalization.
2017
The Gated Recurrent Unit Cho et al. (2014) aggregates the sequence of conditioned variables (of length q) as follows:
u(n]] = gg (u(n 1aan-11; 1), Z, (x. 1)),
(44)
S(so x. f) = D (u(g))
(45)
where [n]] is the hidden memory of a GRU, initialized at zero. FORU denotes the update equations of a GRU we omit gating functions from the notation and D decoder MLP that maps the final GRU hidden state to the desired output, that is, the value of the solution at the desired spatio- temporal coordinate (ar, t), We used two-layered gated recurrent unit with a hidden vector of size 128, and two-layered MLP with 128 neurons activated by the Swish function for D.
Training loop-To create artificial generalization scenarios during training, we employ spatial sub- sampling. Specifically, during each gradient iteration, we randomly and uniformly mask 25% of X and feed the remaining 75% to the output predictor (System 1). To reduce training time further and improve generalization
unseen locations, set of points for querying the model (i.e. as inputs to viy,). To do so, we maintain probability weight vector W of dimension X'' initialized to one. At each gradient descent step, randomly select
use bootstrapping by randomly sampling smaller
16

---

Published as a conference paper at ICLR 2024
N=1, 024 points from X' T weighted by W. We update the weight matrix by setting the values at the sampled locations to zero and then adding the loss function value to the entire vector. This procedure serves two purposes: (a) it keeps track of poorly performing points (with higher loss) and (b) it increases the sampling probability for points that have been infrequently selected in previous steps.
The choice of A in the dynamics loss equation |llloows us to reduce the complexity of the model. In Table we present results obtained with A 3A* indicating that the output predictor (System 1) predicts the latent state representation three time steps later. Consequently, the number of auto- regressive steps during training decreases from T/ (e.g.. for MeshGraphNet and MAgNet) to T/Λ. In Table 2, we used A = 2A" For a more comprehensive discussion on the effect of A on performances, please refer to Appendix G
Training parameters - To be consistent, we trained our model with the same training setup over all different experiments (i.e. same loss function, and same hyper-parameters) However, for the baseline experiments, we did adapt hyper-parameters and used the ones provided by the original authors when possible (see further below). We used the AdamW optimizer with an initial learning rate of 10-3, Models were trained for 4,500 epochs, with a scheduled learning rate decay multiplied by 0.5 after 2,500; 3,000; 3,500; and 4,000 epochs. Applying gradient clipping to a value of effectively prevented catastrophic spiking during training. The batch size was set to 16.
F BASELINES AND DATASETS DETAILS
F.I BASELINES
The baselines trained with the AdamW optimizer with learning rate set at 10-3 for 10,000 epochs on each dataset. We keep the best-performing parameters on the validation set for evaluation on the test set.
DINo - we used the official implementation and kept the hyper-parameter suggested by the authors for Navier and Shallow Water. For Eagle, used the same typper-parameters as for Shallow Water. The training procedure was left unchanged.
implementation of the model in PyTorch, with 8 layers of GNNs for Navier and Shallow Water, and to 15 for Eagle. Other unchanged. We warmed up the model with single-step auto-regressive training with noise injection (Gaussian noise with a standard deviation of 10-4), as suggested in the original paper, and then fine-tuned the parameters by training on the complete available horizon. Both steps try to minimize the mean squared eтror between the prediction and the ground truth. Edges are computed using De- launay triangulation. During evaluation, we perform cubic interpolation between time steps (linear interpolation gives better results on Eagle) first, then 2D cubic interpolation on space to retrieve the complete mesh
MeshGraphNet- we used
were kept
MAgNet- We used our implementation of the MAgNet[GNN] variant of the model, and fol- lowed the same training procedure as for MeshGrapNet. The parent mesh and the query points extracted from the input data using the same spatial sub-sampling technique as in ours, and the edges are also computed with Delaunay triangulation. During evaluation, we split the query points into chunks of 10 nodes, and compute their representation with all the available measurement points. This reduces the number of interpolated vertices in the input mesh and improves performances at the cost of higher computation time (see figure 4 However, to be fair, this increase in computational complexity introduced by ourselves was not taken into account when we discussed computational complexity in appendix
F.2 DATASET DETAILS
Navier & Shallow Water - Both datasets are derived from the ones used in (Yin et al. (022). We adopted the same experimental setup but generated distinct training, validation, and testing sets. For details the GT simulation pipeline, please see Yin et al. (0022). The Navier dataset comprises 256 training simulations of 40 frames each, with additional two times 64 simulations allocated for validation and testing. Simulations are conducted on uniform grid of 64 by 64 pixels (i.e. ?),
17

---

Published as a conference paper at ICLR 2024
Plyssca Space
Space
EC
Romend
Training
Plasicall Space
Latest Souce
LEEET Spice
Jaling Eracneg
slysicl Space
Playsical Space
Evaluation:
SUAC JC
[indacia Qottsiilaoooyplllaa
Figure 4: MaGNet suffers from drastic shifts in distribution between training and evaluation. The model trained on points from X, which corresponds to a small portion of the domain. We used subsampling trick to artificially generate queries. During evaluation. at every available point in the complete simulation, hence, MaGNet must interpolate the initial condition large number of query points, filling the input of the auto-regressive model with noisy estimates of the IC.
require the prediction
measuring the vorticity of fluid subject to periodic forcing. During training, simulations were cropped to 20 frames. The Shallow Water dataset consists of 64 training simulations, along with 16 simulations in both validation and testing. Sequences of length T = 20 were generated. The non-euclidean sampling grid for this dataset is of dimensions 128 64.
Eagle Eagle is a large-scale fluid dynamics dataset simulating the airflow generated by drone within a 2D room. We extract sequences of length T = 10 from examples within the dataset, limiting the number of points to 3,000 (vertices were duplicated when the number of nodes fell below this threshold).
The spatially domn-sampled versions of these datasets (employed in Table and were obtained through masking. We generate a random binary mask, shared across the training, validation, and test sets, remove specified number of points based the desired scenario. Consequently, the observed locations remain consistent across training, validation, and test sets, except Eagle, where the mesh varies between simulations. For Navier and Shallow Water, the High setup retains 25% of the original grid, the Middle setup retains 10%, and the Low setup retains 5%. In the case of Eagle, the High setup preserves 50% of the original mesh while the Low setup retains only 25%. Temporal down-sampling was also applied by regularly removing fixed number of frames from the sequences, corresponding no down-sampling (1/1 setup), half down- wamnlsanling (1/2) and quarter down-sampling (1/4). During evaluation, the models are tasked with predicting the solution to every location and time instant present in the original simulation.
G
MORE RESULTS
Time continuity illustrated in Figure on the Navier dataset. Our model and the baselines are trained in a very challenging setup, where only part of the information is available. During training, not only does the spatial mesh only contains 25% of the complete simulation grid, but also the time-step is increased to four time its initial value. In this situation, the model needs to represent low-resolution data while being trained on sparse data.
Generalization to unseen future timesteps Beyond time continuity, our model offers some gen- cralization to future timesteps. Table sshows extrapolation results for high/mid/low subsampling of the spatial data on the Navier dataset which outperforms the predictions of competing baselines.
Generalization to unseen grid - In our spatial and temporal interpolation experiments (tables and of the main paper), assumed that the observed mesh remains identical during training and testing. Nevertheless, the ability to adapt to diverse meshes is an important aspect of the task. To evaluate this capability, we trained our model in the spatial extrapolation setup on the Navier dataset.
18

---

Published as a conference paper at ICLR 2024
fo +
fo
En + 12
fc 16
tn 20
Goonal Trnii
Ours
DINO
MAgNes MestGmaillNe
Figure 5: Qualitative results Shallow- Shallow-Water Simulation obtained with our model and the baseline in the challenging 5% setup on the Shallow Water dataset (without temporal sub-sampling). Each model is initialized with a small set sparse observations and needs extrapolate the solution at many unseen positions. Our model outperforms the baselines, which struggle to compute the solution outside the training domain.
1
2000
(rag 300
Figure 6: Time continuity on the Navier dataset during training, models only exposed sparse observation of the trajectories, represented spatially by the dots in the upper left figure and temporally by the crmi-aransparent frames. Our model maintains the temporal consistency of the solution and outperforms the baselines.
We compute the error when exposed to different meshes, potentially with a different sampling rate, and report the results in table]4 Our model demonstrates good generalization skills when confronted with new and unseen grids. We observe that the error on new grids close to the error reported in table in the Ext-X' case, show additionally that the model can generalize even if the observed grid is different. Notably, the model performs well when trained with a medium sampling rate. Despite some performance degradation when the evaluation setup is significantly different compared
19

---

Published as a conference paper at ICLR 2024
Navier
High
Mid
Low
In-X'
2.266
2.017
3.154
DINo
Ext-X
2.317
2.136
6.740
In-X
6.853
3.136
1.378
Interp. MGN
Ext-X
7.632
6.890
15.55
In-X
171.5
31.07
10.02
MAgNet
Ext-A
227.0
57.60
89.20
In-X''
0.3732
0.3563
0.3366
Ours
Ext-A'
0.3766
0.3892
0.6520
Table 3: Time Extrapolation We assessed the performances of our model vs. the baselines in a time-xxtrapolation scenario by forecasting the solution on a horizon two times longer than the training one (i.e. 40 frames). Our model remains performant.
Training
Navier
Shallow
High
Mid
Low
High
Mid
Low
In-X
0.2492
0.7929
4.5165
0.5224
1.5431
4.3447
High
Ext-A'
0.2477
0.7782
4.4038
0.5256
1.5822
4.4963
In-X'
0.4370
0.3230
0.9759
0.8528
1.2908
3.6766
Evaluation
Mid
Ext-X
0.4410
0.3401
0,9496
0.8617
1,2589
3.6043
In-X'
2.2000
0.4039
0.6732
2.4395
1.5634
3,4793
Low
Ext-X
2.2037
0.4216
0.7892
2.3914
1.5313
3.2334
Table 4: Generalization to unseen grid We investigate generalization to previously unseen grids by training our model on the Navier dataset in the space extrapolation setup. We report the error (MSE (x 3)) inside and outside the spatial domain X' measured with different sampling rates unseen during training. The diagonal shows results on grids with identical sampling rates wrt. training. but sampled differently. Our model shows great generalization properties.
training, model effectively maintains its interpolation quality between out-of-domain error (Ext-2t) and in-domain error, testifying to the robustness of our dynamic interpolation module.
Ablations study the impact of key design choices in Figure 7p. First, show the effect of the subsampling strategy to favor learning of spatial generalization, c.f. Section 3.4| where we sub- sample the input the auto-regressive backbone by keeping 75% of the mesh. We ablate this feature by training the model on 100%, 50%, and 25% of the input points. When the model is trained on 100% of the mesh, it fails to generalize to unseen locations, as the model is always queried on points lying in the input mesh. However, reducing the number of input points significantly further from the operating point decreases the performance of the backbone, as it does not dispose of enough points to learn meaningful information for prediction. We also replace the final GRU with simpler aggregation techniques, such as a mean and a maximum pooling, which drastically degrades the results. Finally, As expected, this deteriorates the results
ablate the dynamics part of the training loss (Eq. 13 significantly.
More ablation on the interpolator We conducted ablation study show that limiting attention detrimental. To do so, we designed four variants of our interpolation module:
Single attention (w/o GRU) performs the attention between the query and the embed- dings in single shot, rather than time-step per time-step. This variant neglects the insights from control theory presented in section 3.1 (Step 2). The single softmax function limits the attention to a handful points, whereas our method encourages the model to attend to at least one point per time step and reason larger timescale. considering past and future predictions, which is beneficial for interpolation tasks, as supported by proposition 2
Spatial (w/ GRU) & Temporal (w/o GRU) neighborhood - limit the attention to the nearest temporal spatial points, which significantly degrades the metrics. To handle setups with sparse and subsampled trajectories, the interpolation module greatly benefits from not only distant points but also from the temporal flow of the simulation.
20

---

Published as a conference paper at ICLR 2024
Single
Temporal
Spatial
Temporal
ANP Kim et al. (2018)
Ours
attention
attention
neigh.
neigh.
In-X / In-7
0.2113
0.3863
0.2912
0.5623
0.4130
1.734
Ext-X In-T
0.2251
0.4168
0.3180
0.6328
0.6681
1.835
In-X / Ext-T
0.2235
0.4094
0.3095
0.6030
1.9624
1.820
Ext-X / /xx-T
0.2371
0.4388
0.3350
0.6741
2.1818
1.920
Table 5: Ablation on interpolation We performed four ablations on the interpolation module (MSE (x 1x10-3)) Single attention combines all za[ri] into a single key vector, employing attention only once (w/o GRU). Temporal attention replaces the GRU with 2-head attention, Spatial neigh. restricts attention to the five spatially nearest points from the query, and Temporal neigh. computes attention only with the nearest time [r|rz] to the queried time 7 (w/o GRU). These results indicate that considering long-range spatial and temporal interactions is beneficial for the interpolation task.
(a)
Misian Staly
10
Runtime
Mout Squand (ome
(b)
Ao query locaiioss
10F
Runimm
20-1
10
13696
25%
Mean
Mox
Oura
50
100
125
Spocii
Agiegation
#ot query time stagl
Figure 7: Ablations and runtime (a) Ablations Navier (Yn al. 2022 Stokes 2009) with
10% of data and half temporal resolution, from left to right: exploring subsampling strategies, re- placing GRU par mean/max pooling, removing physics grounding. (b) Runtime analysis as a func- tion of query locations and time steps, respectively. The graph shows the average runtime (over 100 runs), and shadows indicate lower and upper bounds over the runs.
Temporal attention (w/o GRU) replaces the GRU in our model with 2-head attention layer. This variant of our model does not improve the performance compared to GRU. We argue that GRU is more suited for accumulating observations time, as its structure matches classic observer designs in control theory.
Attentive Neural Process Kim et al. (2018) is interpolation module close to ours re- sembling the Single attention ablation, with an additional global latent e to account for uncertainties. The model involves prior function g(c, s) trained to minimize the Kullback-
Leibler divergence between @ s sX, (X,T) T)) (computed using the physical state at observed
points) and q( s s X, 10, (0,) (computed using the physical state at query points).
Results are shown in table All ablations exhibit worse performance than ours. Note that the ANP ablation involves performing the interpolation in the physical space to compute the Kullback-Leibler divergence during training. Thus, the interpolation module cannot use the latent space from the auto- regressive module, which may explain the drop in performance, Adaptating ANP directly leverage the latent states is probably possible, but not straightforward and requires several key changes in the architecture.
Efficiency - the design choices we made led to a computationally efficient model, compared to prior work. For all three baselines, the required number of computed time steps for the auto-regressive rollout depends on (1) the number of predicted time steps, and (2) the time values themselves, as for later values of f, more iterations need to be computed. In contrast, our method forecasts using attention from a set of "anchor states", which is controlled through the hyper-parameter A. The length of the auto-regressive rollout is therefore constant and does not depend on the number of predicted time steps. Furthermore, while DINo scales very well to predict additional locations, requires a costly optimization step to compute ag. MGN does benefit from the efficient cubic interpolation algorithm, which is side effect of the fact that it has been adapted to this task, but not
21

---

Published as a conference paper at ICLR 2024
Uza(m
(a) Frontier tracking: when queried on streamline between areas opposite vorticity, the interpolation module attends not only to the spatial neighborhood but also to the temporal flow near the frontier.
8A(x,7 Uza
LURE
(b) Blob tracking: in homogeneous areas, the model tracks the origin of the perturbation, and focuses on its displacement. Our dynamic interpolation exploits the evolution of the state rather than merely averaging neighboring nodes.
LAGE
0 + (2.
EIEL
(c) Periodic boundaries: model effectively leverages the periodic condition of the Navier dataset, espe- cially when queried on points originating from perturbations the other side of the simulation. Again, the interpolation depends on which points explain the output, rather than the neighborhood.
Figure 8: Norm of output derivative wrt. each (x:))(<:)) (Navier, high spatial subsampling setup). We display top-100 nodes (+) with the highest norm, i.c. the most important nodes for interpolation at query point Using gradients rather than attention allows visualize the action of the GRU. We observe context-adaptive behaviors, leveraging temporal flow information over local neighbors, challenging implement in handcrafted algorithms.
designed for it. We experimentally confirm these claims in Figure where we provide the evolution of runtime as a function query locations, and of query time steps, respectively. In both cases, our model compares very favorably to competing methods.
Attention maps - To further support our claims, we analyzed the behavior of the interpolation module in more depth and showed the top-100 most important nodes from the embedding points =a[n](x;) used to interpolate at different queries. The figure is shown Figure S We observed very complex behaviors that dynamically adapt to the global situation around the queried points. Our interpolation module appears to give more importance to the flow rather than merely averaging the neighboring nodes, thus relying on "why" the queried point is in a specific state. Such behavior would be extremely difficult to implement in handcrafted algorithm.
Parameter Sensitivity Analysis We investigate the influence of two principal hyper-parametery namely the step size A and the number of residual GNN layers L. on the performances of our model. We present the results of our experiments figure on the Navier dataset, which has been spatially down-sampled at 10% during training and has a temporal resolution reduced by two.
The choice of the step size between iterations of the auto-regressive backbone directly affects both training and inference time. For a trajectory of T frames, the number of anchor states zu [r] is determined by [T/ Increasing the step size of the learned dynamics leads to higher num- ber of embeddings over which the models need to reason. A parallel can be drawn between this phenomenon and the influence of the discretization size on the accuracy of numerical methods for solving PDEs. Furthermore, the selection of A also impacts the generalization capabilities of the model in Ext-T. When A A the model is queried during training on intermediate instants not directly associated with any of the anchor states za In]. This is visible in Figure9 for instance, with A = A the In-X//nnT error is the lowest, but other metrics increases compared to A = 2
The number of layers L in the rego-regressive backbone significantly influences the overall perfor- mance of the model, both within the domain and on the exteriors. Increasing the number layers generally leads to improved performance. However, it appears that beyond T. the error starts
22

---

Published as a conference paper at ICLR 2024
Infiteence of
Infueence of backhooue depti
1x 10
4 10
3 10--
Squared lieoot
3 10
Squared
2 10
10-1
10"4
10
of lagers
Wedgltting Lihesess
4x
Squard
3.x 10
2 10
10-1
0.001
0.01
0.1
100
Ma lequasiion 13)
(b)
Figure 9: Impact hyper-parameters on model performance We evaluate the impact of three critical hyper-parameters on our architecture, namely, (a) the step size A, the depth L of the auto- regressive backbone (b) the weighting the dynamics cost equation [3] To assess the perfor- mance, employ the 10% Navier dataset with 1/2 frames and compute metrics for both in-domain and out-domain. The results reveal that increasing the depth the GNN layers enhances the model's performance, while lower values of A lead to better metrics. However, we observed degradation in the ability the model to generalize to unseen time instants for the special case A=A". Moreover, found that equally weighting of both terms Leontimaous and Ldynamics leads to best results.
increase, indicating saturation point in terms of performance gain. The relationship between the number layers and model performance is visually depicted Figure Throughout this paper, we maintain this hyper-paramess constant for the sake of simplicity, as our primary focus is the spatial and temporal generalization the solution.
Failure cases we expose failure cases the Eagle dataset (in the Low spatial down-sampling scenario) in figure TO, some particularly challenging instances of this turbulent dataset, we noticed drops in accuracy located in fast-evolving regions of the simulation, and in particular near the flow source. We hypothesize that the origin of the failure is related to the comparatively smaller processor unit used in our auto-regressive backbone compared to the baseline introduced in Janny al. 2023) hence producing less accurate anchor states when the horizon increases.
23

---

Published as a conference paper at ICLR 2024
Lo + 2
Lp + 4
Lp + 6
Lg +
Lo 10
Grocaa
25
20
15
Figure 10: Failure cases on Eagle- We observed failure cases on highly challenging instances of the EAGLE dataset the prediction horizon is increasing. We show the per point error in three different instances and observed that the error increases with the time horizon, especially close to turbulent areas, such as below the UAV.
24


Reviews:
[{"review_id": "IjQou1gaFc", "reviewer": "Reviewer_uqXM", "summary": "In this paper, the authors presented a neural-network-based model that predicts the evolution of a physical system, starting from sparse observations. In order to do so, they consider two (Markovian) dynamical systems (one discrete and one continuous) to predict the time evolution and query the solution in unseen points in the domain.", "questions": "- I would consider adding a bit more information about the weaknesses highlighted before, while they do not change the results, they can provide a more clear picture of the limitations of the model\n- In the introduction, the model from Yin et al. (2022) does satisfy R3, so it would be better to correct the statement there\n- The loss function (eq. 13) has two terms, did you check how the balance of the two terms affects the solution?\n- The author state that there is a trade-off in the performance of $\\phi_q$. Did you train different models to quantify it? What is the dependence of the performance on the training sequence division (q) and length (T). How is it related with the timescales of the physical phenomena tested?\n\n- minor: typo in fig.1, where $\\phi_q$ is $\\phi_p$\n\nAs a final observation, the model performs well but in my opinion the tasks considered are a bit simplistic and do not really highlight the limitations of the proposed models.\nIt would be better to focus on more challenging benchmarks (for instance Biferale et al. https://arxiv.org/abs/2006.07469) in order to bring these models closer to actual applications.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The methodology is well presented and the claims are supported by analytical proofs. The proposed model does provide better results then the models available in the literature, both in time and space", "weaknesses": "While the method works well in practice for the presented datasets, I think that the authors could highlight better the hypothesis behind their model. On the one hand, they correctly highlight how DINo (Yin et al. 2022) assumes the existence of a learnable ODE that models the system dynamics. On the other hand, the authors' approach hypothesizes the existence of a Markovian discrete system that is able to represent the advance in time. This is an equivalently strong hypothesis in my opinion. Furthermore, we assume that an latent space $\\mathbf{z}_d$ can be found using an encoder. This is even more relevant since they perform the time advancement in the latent space: with a scheme \"encode-process-decode\" at each timestep, the encoder needs to take care only of the spatial representation, instead in this case the temporal dynamics also need to be embedded.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0WCVy2ffGP", "reviewer": "Reviewer_eMKs", "summary": "This paper proposes a way to perform forward modeling given sparse initial measurements and is able to retrieve dense solutions out of it. To do so, this paper proposes two systems. One system takes sparse measurement as input - encodes into latent space and propagates in the latent space. Another system takes propagating latent space features and predicts the dense solution. It shows that it works better under sparse observation situations.", "questions": "1. in the webpage - 2D Navier dataset visualization - what does it mean by \"supervision\". Does it mean that the supervision used during the training is those sparse observations in space over time? If the training data is that sparse, how does a data-driven approach learn high-resolution details? \n\n2. as discussed in the paper, sparse observation yields an ill-posed state except the observed time period is long enough. How would the author resolve the ill-posedness? Are there any failure modes observed during the test time? \n\n3. would it be possible to incorporate initial observation over a longer time step (instead of just s_d[0], but maybe s_d[0],s_d[1],...,s_d[n]) to make the problem less ill-posed. Given such a sparse initial state, there could be many plausible initial states, I am curious why would this work...", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Solving an interesting problem where there are only sparse observations, which is common in real-world scenarios. The overall formulation seems novel and there are some interesting theoretical derivations.", "weaknesses": "It is not clear how the intrinsic ill-posedness of the system is addressed - since there is only sparse observation given in the initial state, if i understand correctly.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rvFBuBJnvj", "reviewer": "Reviewer_Kpnn", "summary": "This paper presents a data-driven framework that leverages an auto-regressive structure combined with transformer inference, GNN, and GRU techniques to predict continuous spatial and temporal dynamics in fluid systems from sparse observations. The proposed model demonstrates the ability of generalization to unseen initial conditions.", "questions": "- For System 2, the use of GRU implies iterative discrete predictions rather than time-continuous ones. Does this allow for predictions at arbitrary time frames?\n- If the sampling rate $\\Delta$ does not align with $\\Delta^*$ such that $\\Delta$ is not an integer multiple of $\\Delta^*$, would this lead to enhanced model generalization for unseen time intervals?\n- I am curious about the performance over long-time rollouts. Does executing rollouts in the hidden space contribute to more stable long-term predictions?\n- When the query position $(x, t)$ has a time $t$ that falls between two anchor states, how does the model do inference?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- The paper is well-motivated and well structured\n- The proposed approach is novel to the best of my knowledge, particularly in dividing the dynamics into different systems\n- Including theoretical analysis pertaining to error bounds on two systems enhances the paper's soundness - this approach not only solidifies the work but also offers valuable insights \n- Experiment visualizations are compelling (with nice work on the website), highlighting the advantages of continuous prediction results against relevant baselines in scientific computing", "weaknesses": "- The model's reliance on a regular sampling rate could constrain its applicability. In many practical scenarios, we might not have access to a dataset $\\mathcal{D}$ with a consistent sampling rate. This could impact model training, given that the anchor states are determined based on the available rate.\n- A discussion on how the anchor state computation rate $\\Delta$ influences performance would be beneficial. Given that fluid dynamic systems can exhibit high-frequency changes, a delayed computation rate might fail to capture these rapid transitions.\n\n\n---\n\nThis does not influence my score: what is the meaning of the small Figure on page 9 on the right? Is it a mistake?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "RHRJSWZqTa", "reviewer": "Reviewer_fAqF", "summary": "The paper proposes a data-driven and spatiotemporal continuous surrogate model that generalizes to new unseen initial conditions of partial differential equations. The key idea of the paper is the use of double observable spaces, in which continuous response can be computed on top of learned latent dynamical systems. Training objective is defined on both of the observable spaces. The proposed method is compared with strong baselines in representative experiments.", "questions": "See the weakness above.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is appropriately placed in the current growing literature on scientific machine learning to create a novel spatiotemporal continuous model for tackling complex and wide range of problems in the realm of PDEs. The exposition for the motivation is written crisply and is generally easy to follow. The proposed method is also analyzed from the theoretical aspect and the experiments strengthen the validity of the proposed method.", "weaknesses": "What the authors mean by “Generalize to new unseen initial conditions (ICs)” in R2 is a it unclear. Does this mean initial conditions are even ouf of distribution?\n \n\nIt’s still a bit difficult to assess the degree of the paper’s contribution since the authors miss some literatures relevant to the method proposed in Step1, which are so called “latent unroll” [1, 2]. Although those papers do not provide theoretical results, the idea is essentially same (except the format of input) and its effectiveness is shown empirically. The authors might want to refer those papers in the right place and clarify the contribution of the paper. \n \n\n[1] \"Learning to Accelerate Partial Differential Equations via Latent Global Evolution,\" T. Wu, T. Maruyama, J. Leskovec. NeurIPS, 2022.\\\n[2] “Learning latent field dynamics of PDEs,” A. Sanchez, D. Kochkov, J. A. Smith, M. Brenner, P. Battaglia, and T. J. Pfaff,  NeurIPS, 2020.\n\nThe spacial resolution of the encoder's input in the most experiences (and the number of subsequent samples) is low. It is concerning if the proposed method is scalable to large resolution since $f_1$ is implemented as MPNN that hardly captures global information unless the number of nodes is small. Do you observe any relationship among the number of processor blocks of $f_1$, the spacial resolution, and the resulting error?\n \n\n\n\n**Minor comments**:\n\n* Why does the variance of the runtime of proposed method in Figure 7(b) sporadically get high?\n\n* Typo: \"r1, r2 et r3\" in the last paragraph of the section 1.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bt8zAAcbj1", "reviewer": "Reviewer_VJPj", "summary": "The paper addresses the problem of limitations in traditional physical simulations, which rely on numerical schemes and mesh-refinement methods to balance precision and complexity. The proposed methodology is based on a double observation challenge, which involves extracting regularities from weakly informative physical variables that are sparsely measured in space and time. The authors propose a solution with two interlinked dynamical systems defined on, respectively, the sparse positions and the continuous domain. System 1 is a discrete dynamical model used to compute a sequence of latent anchor states z_d auto-regressively, while System 2 is used to design a state estimator q retrieving the dense physical state at arbitrary locations (x, t). The proposed methodology has several potential implications and effectiveness. First, it allows for predictions in a continuous spatial and temporal domain, which is not possible with traditional methods. Second, it is trained on sparse observations, making it highly adaptable and efficient. Third, the proposed method is generalizable and can be applied to a wide range of physical systems.", "questions": "1. The authors do not provide a thorough comparison of their approach to other data-driven methods in the literature. Can you provide a more detailed comparison of your approach to other data-driven methods in the literature, and how your approach differs from these methods?\n\n2. The paper could benefit from a more detailed analysis of the computational requirements of the proposed methodology. Can you provide more information on the computational requirements of your approach, and how it compares to traditional physical simulations in terms of computational efficiency?\n\n3. The authors do not provide a detailed analysis of the interpretability of the proposed model. Can you provide more information on the interpretability of your approach, and how you plan to address the challenge of interpretability in data-driven methods?\n\n4. The paper could benefit from a more detailed discussion of the potential implications and applications of the proposed methodology beyond the specific examples presented in the paper. Can you provide more information on the potential implications and applications of your approach in other physical systems beyond fluid dynamics?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The proposed methodology represents a significant departure from traditional physical simulations and has the potential to improve the efficiency and accuracy of physical simulations.\n2. The authors demonstrate the effectiveness of the proposed methodology through experiments on three datasets in fluid dynamics and compare their approach to strong baselines, which are outperformed both in classical settings and in the extended new task requiring continuous predictions.\n3. The proposed method is generalizable and can be applied to a wide range of physical systems.\n4. The authors provide a website with an online interactive tool to visualize the results.\n5. The paper is well-organized and clearly presents the proposed methodology and experimental results.", "weaknesses": "1. The paper outlines the utilization of a dynamical system for auto-regressive predictions in a mesh-structured latent space, starting from sparse initial conditions. However, it lacks a comprehensive description of the neural network architecture employed. Clarification is needed on how the model systematically handles and compensates for missing data, as this is a crucial aspect for the practical application and reliability of the proposed method.\n\n2. While the document mentions the use of a specific optimization algorithm for training the model, it omits essential details like the learning rate schedule and batch size. A more thorough examination and justification of the chosen hyperparameters, including regularization parameters and the number of hidden layers, is necessary to ensure replicability and to provide insights into the model's design choices.\n\n3. The paper criticizes traditional physical simulations based on partial differential equations (PDEs) and numerical algorithms for their limitations in complex phenomena and handling sparse observations. It is imperative to acknowledge the extensive success and development of these traditional methods across various applications over the years. Data-driven methods, despite showing potential in overcoming some of these limitations, bring their own set of challenges and limitations. These include the substantial data requirements for training and the often opaque nature of these models, which can make interpretation and validation of results challenging. It would be beneficial for the authors to provide a balanced view, comparing the strengths and weaknesses of both approaches to give readers a comprehensive understanding of the context in which their work sits.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)