Learning the Dynamics of Physical Systems from Sparse Observations with Finite Element Networks
Authors: Marten Lienen
Stephan Günnemann
Venue: ICLR 2022 Spotlight (2022)
Keywords: spatio-temporal
finite
elements
forecasting
continuous
partial
differential
equation
PDE
graph
gnn
time-series

Abstract: We propose a new method for spatio-temporal forecasting on arbitrarily distributed points. Assuming that the observed system follows an unknown partial differential equation, we derive a continuous-time model for the dynamics of the data via the finite element method. The resulting graph neural network estimates the instantaneous effects of the unknown dynamics on each cell in a meshing of the spatial domain. Our model can incorporate prior knowledge via assumptions on the form of the unknown PDE, which induce a structural bias towards learning specific processes. Through this mechanism, we derive a transport variant of our model from the convection equation and show that it improves the transfer performance to higher-resolution meshes on sea surface temperature and gas flow forecasting against baseline models representing a selection of spatio-temporal forecasting methods. A qualitative analysis shows that our model disentangles the data dynamics into their constituent parts, which makes it uniquely interpretable.

Published as a conference paper at ICLR 2022
LEARNING THE DYNAMICS OF PHYSICAL SYSTEMS FROM SPARSE OBSERVATIONS WITH FINITE ELEMENT NETWORKS
Marten Lienen & Stephan Günnemann
Department of Informatics & Munich Data Science Institute Technical University of Munich, Germany (marten.. lienen, guennemann @ in. tum. de
ABSTRACT
We propose new method for spatio-temporal forecasting on arbitrarily distributed points. Assuming that the observed system follows an unknown partial differential equation, we derive a continuous-time model for the dynamics of the data via the finite element method. The resulting graph neural network estimates the instantaneous effects of the unknown dynamics each cell in meshing of the spatial domain. Our model can incorporate prior knowledge via assumptions on the form of the unknown PDE, which induce structural bias towards learning specific processes. Through this mechanism, we derive a transport variant model from the convection equation and show that it improves the transfer performance to higher-resolution meshes on sea surface temperature and gas flow forecasting against baseline models representing a selection of spatio-temporal forecasting methods. A qualitative analysis shows that our model disentangles the data dynamics into their constituent parts, which makes it uniquely interpretable.
INTRODUCTION
The laws driving the physical world are often best described by partial differential equations (PDEs) that relate how a magnitude of interest changes in time with its change in space. They describe how the atmosphere and oceans circulate and interact, how structures deform under load and how electromag- netic waves propagate (Courant & Hilbert, 2008). Knowledge of these equations lets us predict the weather (Coiffier, 2011), build sturdier structures, and communicate wirelessly. Yet, in many cases we only know the PDEs governing a system partially (Isakov, 2006) or not at all, or solving them is too computationally costly to be practical (Ames, 2014).
0
E2
dt
Ø 1
Figure Finite Element Networks predict the instantaneous change of each node by esti- mating the effect the unknown generating dynamics the domain volume that the node shares with its neighbors.
Machine learning researchers try to fill in these gaps with models trained on collected data. For exam- ple, neural networks have been trained for weather forecasts (Shi et al., 2015) and fluid flow simulations (Belbute-Peres et al., 2020), both of which traditionally outcomes of PDE solvers. Even the dynamics of discrete dynamical systems such as traffic (Li et al., 2018) and crowds (Zhang et al., 2017) have been learned from data. A challenge facing these models is the high cost of acquiring training data, so the data is usually only available sparsely distributed in space. Since graphs are natural way structure sparse data, models incorporating graph neural networks (GNNs) have been particularly successful for spatio-temporaa forecasting (Yu et al., 2018; Wu et al., 2019).
Our implementation is available http://www..amm...m.mmddiitttete

---

Published as a conference paper at ICLR 2022
In the domain of physical processes we can reasonably assume that the observed system follows PDE. There are mainly two ways to incorporate this assumption as a-priori knowledge into a model. First, encode a known PDE into a loss function that encourages the model to fulfill the equation (Raissi et al., 2019). Another way to go about this is to derive the model structure itself from known laws such as the convection-diffusion equation (de Bézenac et al., 2018). We will follow the second approach.
Consider dynamical system on a bounded domain S C Id that is governed by the PDE Oeu = (1, x: u, u, u2.,...
on functions u : |0, Q + I... If we have a dense measurement uo Ω2 + R of the current state of the system and solution u that satisfies Eq. (1) for all t E [0, T] and also fulfills the initial condition u(0, 2² = uo (a)) at all points æ E S2, we can use u as a forecast for the state of the system until time T. From a spatio-temporal forecasting perspective, this means that we can forecast the evolution of the system if we have a continuous measurement of the state, know the dynamics F, and can find solutions of Eq. (1) efficiently. Unfortunately, in practice we only have finite number of measurements at arbitrary points and only know the dynamics partially or not at all.
Contributions. An established numerical method for forecasts in systems with fully specified dynamics is the finite element method (FEM) (Brenner et al., 2008). In this paper, we introduce the first graph-based model for opatio-temporal forecasting that is derived from FEM in a principled way. Our derivation establishes a direct connection between the form of the unknown dynamics and the structure of the model. Through this connection our model can incorporate prior knowledge on the governing physical processes via assumptions on the form of the underlying dynamics. We employ this mechanism to derive specialized model for transport problems from the convection equation. The way that the model structure arises from the underlying equation makes our models uniquely interpretable. We show that our transport model disentangles convection and the remainder of the learned dynamics such as source/sink behavior, and that the activations of the model correspond to learned flow field, which can be visualized and analyzed. In experiments on multi-step forecasting of sea surface temperature and gas flow, our models are competitive against baselines from recurrent, conporal-oonvolutional and continuous-time model classes and the transport variant improves upon them in the transfer to higher- resouutin meshes.
2
BACKGROUND
2.1 FINITE ELEMENT METHOD
In the following, we will outline how approximate solution the dynamics in Eq. (1) from an initial value ug by discretizing u in space using finite elements. Let X' be a set of points with a triangulation T of d-dimensional, non-overlapping simplices
X = (x(1) E RdM
A C 22 1 1d+
(2)
such that USET CH(A) A) equals the domain l where CH(A) is the convex hull of simplex A. So we define simplex A E T representing the j-th mesh cell as the set of vertices of the cell and denote the domain volume covered by the cell by the convex hull CH (i)) of the vertices. We will assume be scalar field, i.e. u : [0, l IR. If u is vector field, treat it as a system of m scalar fields instead. For a detailed introduction to FEM, we refer the reader to Igel (2017).
Basis Functions. A priori, we assume that the unknown solution u to our problem lies in an infinite- dimensional function space 24 The first step in FEM to make the problem numerically feasible is to approximate 24 with dinee-dameniionll linear subspace L. This subspace can then be written in terms of linear combinations of basis functions i = span {1(1) ... (N))) There are many possible bases and the choice determines various qualities of the resulting procedure such as continuity of the approximation and the sparsity pattern of the mass matrix in Eq. (7).
In our case, we choose the so-called P1 basis of piecewise linear functions (hat functions), see Fig. 2a (Igel, 2017). There are as many basis functions as there are points and each is uniquely defined by being linear when restricted to single cell A T and the constraint
[1
if (6) (0)
(2())) = ( otherwise
(z(() X.
(3)
2

---

Published as a conference paper at ICLR 2022
Initial state ug Approximation
Term .pp (8)
0
0
Timo
21-2
34-1
Ti+1
Ji+2
31
322
273
J74
375
36
277
21
3/2 23
3.4
35
76
(a) Choice basis functions
(b) Discretization state and dynamics
(c) Solving forward in time
Figure 2: Solving a PDE with the Galerkin method and method of lines consists of three steps.
So the basis function 600) is 1 at σ6 falls linearly on mesh cells adjacent to x(1) and is everywhere else. The resulting finite-dimensional function space úf is the space of linear interpolators between values at the vertices, see Fig. 2b. An important property is that if we expand d in this basis, the value of u at the i-th node is just its i-th coefficient.
u(x(())= =
-E, j=1 (x))) =
Galerkin Method. A piecewise linear approximation E ú is not differentiable everywhere and therefore cannot fulfill Eq. (1) exactly. So instead of requiring an exact solution, we ask that the residual R(u)= = 0 u - F (t, x, u, ...) be orthogonal to the approximation space i with respect to the inner product (u. 2) Ja uxx) v(())) dæ at any fixed time t. In effect are looking for the best possible solution U. Because decomposes into N equations, one for each basis function.
generated by finite basis, the orthogonality requirement
(R(u), v/a = 0 Ve EU
(R(u), (l) )n 0
Vi I,.., N
(5)
Plugging the residual back in and using the linearity of the inner product, we reconstruct system of equations that resemble the PDE that we started with.
(0.u, ( = (F (t, u....,,)0))))) a
vi=
(6)
At this point we can stack the system of N equations into vector equation. If we plug in the basis expansion Σド=1 (5) for into the left hand side, get linear system
A0c-= m
where Aij = (0p(0) (000),00))))) is the called mass matrix, e is the vector of basis coefficients of u, and my (F' (t, u.... (p()) captures the effect of the dynamics F. The left hand side evaluates to Ad,c, because the basis functions are constant with respect to time. The right hand side cannot be further simplified without additional assumptions on
Method of Lines. If we can evaluate the right hand side m, we can solve the linear system in Eq. (7) for the temporal derivatives of the coefficients of at each point time. fact we have converted the PDE into a system of ordinary differential equations (ODEs) which we can solve with an arbitrary ODE solver given an initial value ( as in Fig. 2c. This is known as the method of lines because we solve for along parallel lines in time.
To find a vector field u [0, 7] x Ω > RT instead of a scalar field, we treat the m dimensions of u as system of scalar fields. This results m copies Eq. (7), which we need to solve simultaneously. Because the mass matrix A is constant with respect to we can combine the system into a matrix equation
- -
where C,M
M E Nxm are the stacked e and m vectors, respectively. In summary, the spatial discretization with finite elements allows us to turn the PDE (1) into the matrix ODE (8).
2.2 MESSAGE PASSING NEURAL NETWORKS
Pesssage-Passingg Neural Networks (MPNNs) are a general framework for learning on graphs that encompass many variants of graph neural networks (Gilmer al., 2017). It prescribes that nodes in
3

---

Published as a conference paper at ICLR 2022
graph iteratively exchange messages and update their state based on the received messages for P steps. For a graph g (V.E) with nodes V and edges E, and initial node states h Ve V, the p-th propagation step is
h Supd hy (11-1), E (U, Smsg (n.e (h n nn 1)
where Smsg maps node states and edge attributes to messages and fupd updates node's state with the aggregated incoming messages. The final node states h (P) can then be interpreted as per-node predictions directly or passed as node embeddings downstream systems.
In this work, we employ a slight generalization of the above to undirected hypergraphs, i.c. graphs where the edges are sets of an arbitrary number of nodes instead of having a cardinality of exactly For such hypergraph G = (V,E) E) with nodes V and hyperedges {u, w... E E, and initial node states h. Ve E V. the p-th propagation step is
h (p) fupa h
(11-11),
Σ
fmag (
(10)
CEE s.t.VEe
Note that Smsg jointly computes a separate message for each node v participating in a hyperedge e.
3 FINITE ELEMENT NETWORKS
N Consider set nodes X' = (xx(1) (R= node we measure m features y (fo..) E Thm" Our goal is to predict ( at future timesteps i, by solving PDE (1). For it to be well-defined, we need continuous domain encompassing the nodes X.
Ji=l representing N points in d-dimensional space. At each
By default, we construct such a domain by Delaunay triangulation. The Delaunay algorithm triangulates the convex hull of the nodes into set of mesh cells d-dimensional simplices
T = ( C X || d + 1 which we represent as sets of nodes denoting the cell's ver- tices. The PDE will then be defined on the domain f2 U3e7 CH(A). One shortcoming of this algorithm is that it can produce highly acute, sliver-like cells on the boundaries of the domain if interior nodes are very close to the boundary of the convex hull. We remove these slivers post-processing step, because they have small area and thus contribute negligible amounts to the domain volume, but can connect nodes which are far apart. See Appendix H for our algorithm. As alternative Delaunay triangulation, mesh can be explicitly specified to convey complex geometries such as disconnected domains, holes in the domain and non-convex shapes in general.
As first step towards solving the PDE, we need to represent y (fo) as a function on Ω. We define the function ux(x); = Σ Vi (ta) (z)) such that its coefficients as vector in the P1 basis encode the data and evaluating it at the nodes uo (+) (fo0.) agrees with the data by construction.
In Section 2.1, we have seen that solving PDE (1) approximately under an initial condition u(lo, .r) = (x( can be cast as an initial value ODE problem over the trajectory of the coefficients of solution in time. Above, we encoded the features (fo) in the coefficients of Mg. meaning that Eq. (8) exactly describes the trajectory the features through time. The dynamics of the features are thus given by
AO, Y (0) M
(11)
where (())
in yk is the feature matrix and Aij (1)()) so (o()) )) is the mass matrix. We call the matrix M on the right hand side the message matrix and it is given by
(t.5)
Ma = (F (t, x. U, uy...) s()))
(12)
where u encodes the predicted features Y at time in its coefficients.
Our goal is to solve the feature dynamics (11) with an ODE solver forward in time to compute the trajectory of Y and thereby predict (tt)), To apply ODE solver, have to compute o,Y at each time which in turn consists of, first, evaluating the right hand side and, second, solving the resulting sparse linear system.
4

---

Yet, if we were to introduce a deep model fa Rs F in Eq. (14), we would run into two problems. First, each evaluation would require the numerical integration of neural network over each mesh cell against multiple basis functions and adaptive-step ODE solvers can require hundreds of evaluations. Second, conditioning fe on just t, and u(x) would not provide the model with any spatial information. Such a model would not be able to estimate spatial derivatives 0r internally and could Published as a conference paper at ICLR 2022
Solving Eq. (11) for a,Y with the exact sparse mass matrix A carries performance penalty GPUs, because sparse matrix solving is difficult to parallelize. To avoid this operation, we lump the mass matrix. Lumping of the mass matrix is standard approximation in FEM for time-dependent PDEs to reduce computational effort (Lapidus & Pinder, 1999), which leads to good performance in practice.
A.
(0 (0
(13)
We employ the direct variant that diagonalizes the mass matrix by row-wise summation, making the inversion of A trivial.
The remaining step in solving Eq. (11) evaluating the right hand side. Let T
{A E T x00) A} be the set of all mesh cells that are adjacent to node 201), Then we can
break the elements of M up into the contributions from individual cells adjacent to each node.
Mas = (F (t, x, U. ...); (p(1)) ) - ESET, (F (f, (t, It u. u....) (0()) (CH())
(14)
The sum terms capture the effect that the dynamics have on the state the system such as convection and chemical reaction processes given the geometry of the mesh cell, the values the nodes, and the basis function to integrate against. However, the dynamics F are of course unknown and need to be learned from data.
therefore only represent PDEs that are actually just ODEs.
We deal with both problems at once by factoring the inner products in Eq. (14) into
(F (t, x. u, ...) Ap JCH(A)))))))))))))))) = F(. (1, sp (CH())
(15)
Such scalar coefficient (()) A. always exists, because (1, ())))) is constant with respect to and and never O. On the basis of Eq. (15), we can now introduce a deep model So Rs F(0) With this, we sidestep the numerical integration of a neural network, because F(O) constant per mesh cell and the coefficients (I, (L,)()))C())))) can be precomputed once per domain. At the same time, 1-(0)) depends jointly all æ in mesh cell as well as their solution values u(x). Therefore our deep model fe naturally also needs to be conditioned on these spatially distributed values, giving it access to spatial information and making it possible learn actual PDE dynamics.
3.1 MODEL
With these theoretical foundations in place, we can formally define Finite Element Networks (FENs). Let = (.Y, T) the undirected hypergraph with nodes X' and hyperedges T corresponding to the meshing of the domain, so each mesh cell A E T forms a hyperedge between its vertices. Let fe be a neural network that estimates the cell-wise dynamics FA from time, cell location, cell shape and the values at the vertices
(t.1) f0.2 = (t.HA TA (1)
F(
(16)
(+)) E A} are the local coordinates of the where μa is the center of cell A, CS = (x(1) us
cell vertices and ز (t)) (t) }} are the features at the vertices time t. We call fo the free-form term, because it does not make any assumptions on the form of F and can in principle model anything such as reactions between features sources and sinks. fo conditions the cell centers μ and the local vertex coordinates øA separately to help the model distinguish between cell position and shape, and improve spatial generalization.
Define the following message and update functions in the MPNN framework
Smsg (A) (t..) (1,5(1) (CH(A) (1,H()))H(()))
fund m)()) A 1
(17)
5

---

Published as a conference paper at ICLR 2022
where m ΣNET fisg (A),..) ai are the aggregated messages for node x(). Performing one message passing step (P= 1) in this model is exactly equivalent to solving for the feature derivatives 0, ( in Eq. (11) with the lumped mass matrix from Eq. (13) and the message matrix in Eq. (14) with the factorization from Eq. (15) and fa as model for Fa. So making a forecast with FEN for T points in time t = (ta, f1...,, ir) based on sparse measurement y means solving an ODE where at each evaluation step l, we run one message passing step the MPNN described above to compute the current time derivatives of the features k-mt
y ((o. = ODESolve(y (to), O.y, t)
(18)
Because FENS model the continuous-time dynamics of the data, they can also be trained on and make predictions irregular timesteps.
By controlling the information available to fe. we can equip the model with inductive biases. If we omit the current time t, the learned dynamics become necessarily autonomous. Similarly, if do not pass the cell position μ the model will learn stationary dynamics.
An orthogonal and at least as powerful mechanism to induce inductive biases into the model is through assumptions on the form of the unknown dynamics F. Let us assume that the dynamics are not completely unknown, but that a domain expert told us that the generating dynamics of the data include convective component. So the dynamics for feature & are of the form
F (1, x U, ... = V (o (t,, æ u. u...)u)) + F' (u, Og
(19)
where (*) 2, ...) Fd is the divergence- freee velocity field of that feature and F" represents the still unknown remainder of the dynamics. We assume the velocity field to be divergence-free, because we want the velocity field to model only convection and absorb sources and sinks in the free-form term.
By following similar steps with the convection term - (u( x, u, ...))) as for the unknown dynamics, we arrive specialized term for modeling convection. Let go be neural network that estimates one velocity vector per cell and attribute.
(t.i) = go (t.HA.RA, = vtt, x, u.... E Rmd 3
(20)
Ensuring that the learned velocity field is globally divergence- would be costly operation, but by sharing the same predicted velocity vector between all vertices of each cell, we can cheaply guarantee that the velocity field is at least locally divergence free within each cell. The corresponding message function, which we derive in Appendix B, is
(
Jmsg (A)
(1: (A)
(t.i) Σ Eavies (1 + (00
(Vibt), (0 )CH(A))
(21)
With this we can extend FEN with specialized mradsport-modeling capabilities by adding the two message functions and define Transport-FEN (T-FEN) as a FEN with the message function
TFEN mmg
(A)a. = fman (A) fmag (A) (1))
(22)
Consequently, T-FEN learns both velocity field to capture convection and a free-form term that fits the unknown remainder of the dynamics.
Network Architecture. We instantiate both fa and go as multilayer perceptrons (MLPs) with tanh non-linearities. The input arguments are concatenated after sorting the cell vertices by the angle of (+) 14 in polar coordinates. This improves generalization because fo and go do not have to learn order invariance of the nodes. Both MLPs have an input dimension of + + (d + 1)- (m + d) in the non-autonomous, non-stationary configuration since | d + 1. The free-form term fo computes one coefficient per vertex and attribute and therefore has an output dimension (d+ 1) mm whereas the transport term estimates a velocity vector for each cell and attribute resulting in an output dimension of m d. The number of hidden layers and other configuration details for each experiment listed in Appendix D. We use scikit- fem (Gustafsson & McBain, 2020) compute the inner products between basis functions over the mesh cells and the dopri5 solver in tor torchdie ffe (Chen et al., 2018) to solve the resulting ODE. See Appendix E for details model training.

---

Published as a conference paper at ICLR 2022
4 EXPERIMENTS
Baselines. We evaluate FEN and T-FEN against selection of models representing variety of temporal prediction mechanisms: Graph WaveNet (GWN) combines temporal and graph convolutions Difference Graph Network (PA-DGN) estimates spatial derivatives
(Wu et al., 2019); Physics- as additional features for a recurrent graph network (Seo et al., 2020); the continuous time MPNN- based model by (lakovl et al., 2021) which we will call CT-MPNN, uses general MPNN to learn the continuous-time dynamics of the data.' See Appendix D for the configuration details.
GWN and PA-DGN use 12 and 5 timesteps respectively as input while the ODE-based models make their predictions based on a single timestep. To ensure a fair comparison, we do not evaluate any of the models on the first 12 timesteps of the test sets.
Datasets, For our experiments we have chosen two real-world datasets and a synthetic one. The Black Sea dataset provides data on daily mean sea surface temperature and water flow velocities on the Black Sea over several years. ScalarFlow is collection of 104 reconstructions of smoke plumes from seconds long, multi-view camera recordings (Eckert et al., 2019). In each recording fog machine releases fog over heating element in the center of the domain which then rises up with the hot air through convection and buoyancy. CylinderFlow contains simulated fluid flows around a cylinder along channel simulated with the inviscid Navier-Stokes equations (Pfaff et al., 2021).
All our experiments are run on hold-out datasets, i.e. the year 2019 of the Black Sea data and the last 20 recordings of ScalarFlow. See Appendix C for detailed description of the datasets, our sparse subsampling procedure, choice of normalization, and mnin--alidation-tee test splits.
Table 1: Mean absolute error and number of function evaluations for multi-step forecasting on the test set. We forecast 10 steps ahead and averaged the metrics over 10 runs. Results marked with include models that could not be trained to completion because of excessive memory requirements.
ScalarFlow
Black Sea
CylinderFlow
MAE x
NFE
MAE xII
NFE
MAE xIO
NFB
PA-DGN
13.86 + 0.14
9.20 = 0.05
5.46 = 0.02
GWN
9.79 # 0.05
8.94
3.04 = 0.28
CT-MPNN
8.97 # 0.06
157.2
8.54 =.54 0.05
32.9
5.18 = 0,43+
42.2
FEN
9.14
23..0=43.0
8.78 0.11
102.3
2.87 = 0.17
123.6 33.1
T-FEN
9.04 0.06
315.5 = 165.0
6..4
4.09 = 0.62+
67.1 + 22.4
Multi-step Forecasting. The main objective of our models as well as the baselines are accurate forecasts and we train them to minimize the mean absolute error over 10 prediction steps. We chose this time horizon because it is long enough that a model needs to learn non-trivial dynamics to perform well and more accurate dynamics lead to better performance. At the same time it is short enough that not too much uncertainty accumulates. For example, if we train any of the models on Black Sea on much longer training sequences, they learn only diffusive smoothing because the correlation between the temperature distribution today and in 30 days is basically nil, so that the only winning move is not too play - similar to how forecasting the weather for more than 10 days is basically impossible (Mailier, 2010).
For Table we have evaluated the models on all possible subsequences of length 10 of the respective test set of the ScalarFlow and Black Sea datasets and on all consecutive subsequences of length 10 in the CylinderFlow test set because of the size of the dataset. The results show that FEN is competitive with the strongest baselines on all datasets. On real-world data, the transport component in T-FEN provides a further improvement over the base FEN model. The number of function evaluations for all ODE-based models depend strongly on the type of data. While CT-MPNN is the cheapest model on ScalarFlow in terms of function evaluations by the ODE solver, FEN and T-FEN are more efficient on Black Sea. On the synthetic CylinderFlow dataset, neither CT-MPNN nor T-FEN can by trained fully because the learned dynamics quickly become too costly to evaluate. These results show that including prior knowledge on the form of the dynamics in the model can improve the prediction error when the additional assumptions hold.
IGWN and PA-DGN condition their predictions on multiple timesteps, giving them a possible advantage over ODF-hased models since they can the additional input for denoising or to infer latent temporal information.
7

---

Published as a conference paper at ICLR 2022
CT-MPNN
FRN
T-PRN
Data
FEN
T-FEN
CT-MPNN
0.95
Black Sea
0.9
to + 805ms
MAE
ScalurFlow
0.09
3000
5000
Nodes
Figure 3: Predictive accuracy of mod- els trained on 1000 nodes and evalu- ated on increasingly fine subsamplings of the test data.
Figure 4: Evolving the learned dynamics of FEN models and the strongest baseline trained on 10-step forecasting forward for 60 time steps reveals the differences between these models to the human eye.
Cuper-Resolution. Graph-based spatio-temporal forecasting models make predictions on all nodes jointly. Training these models on many points or fine meshes therefore requires a prohibitive amount of memory to store the data points and large backwards graphs for backpropagation. One way around this problem train on coarse meshes and later evaluate the models on the full data. In Fig. 3 we compare FEN and T-FEN against the strongest baseline on the Black Sea dataset. While all three models increase in prediction error as the mesh resolution increases, FEN models deteriorate at a markedly slower rate than CT-MPNN. We ascribe this to the fact that in FEN and T-FEN estimation of the dynamics and its effect on the nodes are to a certain extent separated. The terms fo and ge estimate the dynamics each cell, but the effect the node values is then controlled by the inner products of basis functions which incorporate the shape of the mesh cells in a theoretically sound way. In CT-MPNN on the other hand, the model learns these aspects jointly which makes it more difficult to generalize to different mesh structures. See Appendix F for the complete results on both datasets and visual comparison of low and high resolution data.
Extrapolation. While we train the models on 10-step forecasting, it is instructive to compare their predictions on longer time horizons. Fig. 4 shows the predictions of the three strongest models smoke plume after 60 steps. Because these models are all ODE-based, they were conditioned on single timestep and have been running for hundreds of solver steps this point. First, we notice that all three models have the fog rise at about the same, correct speed. However, with CT-MPNN the fog disappears at the bottom because the model does not learn represent the fog inlet as a fixed, local source as opposed to our models. Comparing FEN and T-FEN, we see that the former's 's dynamics also increase the density further up the fog column where no physical sources exist, while the latter with its separate transport term modeling convection keeps the density stable also over longer periods of time. See Appendix G for a comparison with all baselines at multiple time steps.
Interpretability. In FEN the free-form term fe models the dynamics F' black box and as such its activations are equally opaque as the dynamics we derived from. However, as we have seen in Section 3.1, by making assumptions on the structure of F we impose structure onto the model and this structure in turn assigns meaning to the model's ' components. For a T-FEN we see that the transport term corresponds remainder of the dynamics. Fig. 5 shows that this disentanglement of dynamics works in practice. In the ScalarFlow data, we have a fog inlet at the bottom, i.c. source of new energy, and from there the fog rises up by convection. Plotting the contributions of fe and to 0,y separately shows that the free-form term represents the inlet while the transport term most active in the fog above. Instead of inspecting the contributions to i.y, we can also investigate the learned parameter estimates directly. Fig. 6 shows that the transport term learns a smooth velocity field for the sea surface temperature on the Black Sea.
velocity field for each attribute, while the free-form term absorbs the
5 RELATED WORK
Spatio-Temporal Forecasting, Various DNNs, and in particular GNNs, have been designed to model complex spatial and temporal correlations in data, with applications ranging from forecasting traffic
8

---

Published as a conference paper at ICLR 2022
Free-Form
Transport
Figure 5: Contribution to dyy from the free-form and transport term in T- FEN on a snapshot of ScalarFlow.
Figure 6: Learned temperature flow field in T-FEN after 10 steps on Black Sea data. The model recognized the physical relationship between the features.
flow (see (Jiang & Luo, 2021) and (Yin et al., 2021) for a comprehensive survey), to forecasting COVID-19 (Kapoor al., 2020) and sea-surface temperature (de Bézenac al., 2018). For detailed review of both classical and deep learning approaches, which can be roughly categorized as based on point-cloud trajectories, regular grids, or irregular grids, see Shi Yeung (2018). As one representative of GNN- based methods, which tend to be SOTA, we compare FENs with GWN (Wu al., 2019). GWN uses combination of graph convolution and 1D dilated temporal convolutions to learn spatio-temporal properties of traffic networks.
Neural Networks and Differential Equations. In 1998 Lagaris et al. (1998) introduced the first neural network approach to solve ODEs and PDEs. They formulate the optimization goal as fulfilling the properties the differential equation. Raissi et al. (2019) extend this idea and propose framework to respect any given general nonlinear PDE for continuous- and discrete-time models. Another neural approach to learn PDEs from data is PDE-Net (Long et al., 2018). Its two main components are: learning the differential operators of PDEs through convolutional kernels and applying a neural network to approximate the nonlinear response functions. Building on that work, PDE-Net 2.0 (Long al., 2019) uses symbolic network to approximate the nonlinear response function. de Bézenac et al. (2018) derive a model from the convection -iiffiitiooniiiiiiion equation and estimate the flow field a regular grid with CNN. Chen et al. (2018) leverage ODE solvers their neural architecture to create model family of continuous depth called neural ODEs. Ayed et al. (2019) build upon neural ODEs make spatio-tempon points jointly with the structure of the domain encoded via the node positions as additional features. GNODE (Poli al., 2019) brings the framework neural ODEs of Chen et al. (2018) into the graph domain and provides continuous depth GNNs to make traffic predictions. In dynamical systems the variable depth component coincides with the time dimension which allows the network to work with observations at varying timesteps. Iakovl et al. (2021) extend GNODE with distance vectors between nodes as edge features and use it for PDE prediction. Mitusch et al. (2021) derive a model from FEM similar our work but integrate the learned dynamics against the basis functions numerically and disregard the mesh to avoid designing the discretization into the model.
forecasts by learning the dynamics of all data
6 CONCLUSION
We have introduced Finite Element Networks as new spatio-temporal forecasting model based on the established FEM for the solution of PDEs. In particular, we have shown how FENs can incorporate prior knowledge about the dynamics of the data into their structure through assumptions on the form of the underlying PDE and derived the specialized model T-FEN for data with convection dynamics. Our experiments have shown that FEN models are competitive with the strongest baseline model from three spatio-temporal forecasting model classes on short-term forecasting of sea surface temperature and gas flow and generalize better to higher resolution meshes at test time. The transport component of T-FEN boosts the performance and additionally allows meaningful introspection the model. Its structure is directly derived from the assumed form of the underlying PDE and as such the transport and free-form terms disentangle the dynamics into a convection component and the remainder such inlets and sinks.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
We thank Leon Hetzel for helpful discussions and suggestions.
This research was supported by the Bavarian State Ministry for Science and the Arts within the framework of the Geothermal Alliance Bavaria project.
ETHICS STATEMENT
We propose general model for spatio-temporal forecasting which can be applied in variety of settings. Among these are traffic and crowd flow predictions which have the potential for abuse. However, due to the close connection PDEs, see primarily as a model for physical systems in earth sciences and engineering, with impacts in the study of climate change and understanding of complex systems. Finally, we can also imagine its application in modeling the spread of infectious diseases which also employs PDEs today (Viguerie et al., 2021).
REPRODUCIBILITY
To maximize the reproducibility our experimental results, have fixed seeds for every domain subsampling, training run and evaluation and recorded these in configuration files, so that everything can be re-run. We obtained all seed values from numpy random. SeedSequencc to guarantee that the streams of random numbers from the random number generators (RNGs) are independent and have high sample quality. While we fixed all seeds, model training is still non-eeteriiiistic because the message passing step in graph-based models relies on scatter operation. Scattering implemented via atomic operations on GPUs, which can re-order floating point operations such additions between runs with the same seed inducing non-determinism (Fey & Lenssen, 2019). Nonetheless, we have observed that multiple runs from the same seed only diverge slowly as training progresses and the variability in performance between models anyway as can be seen in the standard deviation of the mean absolute error in Table 1.
from different seeds is small for all
We have used NVIDIA GeForce GTX 1080 Ti for all experiments and evaluations.
REFERENCES
William Ames. Numerical methods for partial differential equations. Academic Press, 2014.
Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari. Learning dynamical systems from partial observations. CoRR, abs/ abs/1902.1 1136, 2019.
Filipe de Avila Belbute-Peres, Thomas Economon, and Zico Kolter. Combining differentiable pde solvers and graph neural networks for fluid flow prediction. In ICML, 2020.
Susanne C Brenner, L Ridgway Scott, and Ridgway Scott. The mathematical theory of finite element methods, volume 3. Springer, 2008.
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In NeurIPS, 2018.
Jean Coiffier. Fundamentals of numerical weather prediction. Cambridge University Press, 2011.
Richard Courant and David Hilbert. Methods of Mathematical Physics: Partial Differential Equations. John Wiley & Sons, 2008
Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical processes: Incorporating prior scientific knowledge. In ICLR, 2018.
Marie-Lena Eckert, Kiwon Um, and Nils Thuerey. Scalarflow: large-scale volumetric data set of real-world scalar transport flows for computer animation and machine learning. ACM Transactions Graphics, 3-:1-16, 2019.
10

---

Published as a conference paper at ICLR 2022
Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural
message passing for quantum chemistry. In ICML, 2017.
Tom Gustafsson and G. D. McBain. scikit-fem: A python package for finite element assembly. Journal of Open Source Software, 5:2369, 2020.
Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre
Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585:357-362, September 2020.
D. Hunter. Matplotlib: A 2d graphics environment. Computing in Science Engineering, 9:90-95, 2007.
Valerii lakovi, Markus Heinonen, and Harri Lähdesmäki. Learning continuous-time pdes from sparse data with graph neural networks. In ICLR, 2021.
Heiner Igel. Computational Seismology: A Practical Introduction. Oxford University Press, 2017.
Victor Isakov. Inverse problems for partial differential equations, volume 127. Springer, 2006.
Weiwei Jiang and Jiayun Luo. Graph neural network for traffic forecasting: A survey. CoRR, abs/2101.11174, 2021.
Amol Kapoor, Xue Ben, Luyang Liu, Bryan Perozzi, Matt Barnes, Martin Blais, and Shawn O' O'Banion. Examining COVID-19 forecasting using spatio-temporal graph neural networks. CoRR, hoss2000.03111. 2020.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Georg Kohl, Kiwon Um, and Nils Thuerey. Learning similarity metrics for numerical simulations. In ICML, 2020.
Isaac Lagaris, Aristidis Likas, and Dimitrios Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. Neural Networks, IEEE Transactions on, pp. 987 1000, 10 1998.
Leon Lapidus and George F Pinder. Numerical solution of partial differential equations in science and engineering. John Wiley & Sons, 1999.
Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. In ICLR, 2018.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In ICML, 2018.
Zichao Long, Yiping Lu, and Bin Dong. Pde-net 2.0: Learning pdes from data with numeric- symbolic hybrid deep network. ArXiv, abs/1812.04426. 2019.
Pascal J Mailier. Can we trust long- range weather forecasts? Management of weather and climate risk the energy industry, pp. 227-239, 2010.
Met Office. Cartopy: a cartographic python library with a matplotlib interface. Exeter, Devon, 2010 2015. URL http: http://ccccoooossoo org. w.l.kk//aaaa hic//rrrpy....
Sebastian K. Mitusch, Simon W. Funke, and Miroslav Kuchta. Hybrid FEM-NN models: Combining artificial neural networks with the finite element method. CoRR, abs/2101.00962,, 2021.
Derek Onken and Lars Ruthotto. Discretize-optimizz vs. optimize- discretize for time-series regression and continuous normalizing flows. CoRR, abs/2005. 13420. 2020.

---

Published as a conference paper at ICLR 2022
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury. Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Des- maison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Te- jani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang. Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In NeurIPS. Curran Associates, Inc., 2019. URL http://papers, neurips cc pee/papr///// 9015 -art--hptrcc-h-c an imperative style -ieeehighh performance deep- seeningg-lbbrrry pdf.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W. Battaglia. Learning mesh- based simulation with graph networks. In ICLR, 2021.
Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. CoRR, abs/1911.07532, 2019.
Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics- innfire- neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Alex Rogozhnikov. Einops: Clear and reliable tensor manipulations with einstein-like notation, In ICLR, 2022.
Sungyong Seo, Chuizheng Meng, and Yan Liu. Physics-aware difference graph networks for sparsely- observed dynamics. ICLR, 2020.
Xingjian Shi and D. Yeung. Machine learning for spatiotemporal sequence forecasting: survey. ArXiv, abs/1808.06e 06865, 2018.
Xingjian Shi, Zhourong Chen, Hao Wang, Din Yeung, Wai-Kin Wong, and Wang-chun Woo. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. In NeurlPS, 2015.
Sam Shleifer, Clara McCreery, and Vamsi Chitters. Incrementally improving graph wavenet perfor- mance traffic prediction. CoRR, abs/1912.07339, 2019.
Alex Viguerie, Guillermo Lorenzo, Ferdinando Auricchio, Davide Baroli, Thomas JR Hughes, Alessia Patton, Alessandro Reali, Thomas E Yankeelov, and Alessandro Veneziani. Simulating the spread of covid-19 via spatially-resolved (seird) model with heterogeneous diffusion. Applied Mathematics Letters, 111:106617, 2021.
Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chenggi Zhang. Graph wavenet for deep spatial-temporal graph modeling. In I/CAI, 2019.
Omry Yadan. Hydra framework for elegantly configuring complex applications. Github, 2019.
URL https:://ttthhiiee wwtp:///ceeccckrsserr///rooeedoe
Xueyan Yin, Genze Wu, Jinze Wei, Yanming Shen, Heng Qi, and Baocai Yin. Deep learning traffic prediction: Methods, analysis and future directions. IEEE Transactions on Intelligent Transportation Systems, 2021.
Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. In IJCAI, 2018.
Junbo Zhang, Yu Zheng, and Dekang Qi. Deep spatio-temporal residual networks for citywide crowd flows prediction. AAAI, 2017.
12

---

Published as a conference paper at ICLR 2022
A DISCUSSION
We have presented a way to derive a spatio-temporal forecasting model directly from FEM. In addition to establishing connection between this well-known theory behind many scientific simulations and neural ODE models, our derivation opens a way to include knowledge in the form of PDEs directly into the structure of machine learning model, as we have shown by the example of the convection equation and T-FEN Extensions of the model with other first or second order, linear PDEs such as the heat equation or the wave equation should follow with similar derivation. However, our approach cannot be applied directly when derivatives of higher order are involved.
The problem fundamentally stems from our choice of basis functions. The P1 basis piecewise linear and thus has vanishing second and higher order derivatives. This means that, at the end of the model derivation, any basis function involved in fixed weights such as (go0), ()))) can appear with at most first derivative, otherwise the expression will evaluate to O. A solution would lie in extending our approach to higher order basis functions, though this is a non-trivial endeavor. In such a setting, the values the nodes of a triangular mesh would no longer suffice to determine all coefficients in the new basis. In future research, we will work on constraining these new coefficients and make learning with higher-order basis functions a well-posed learning problem.
B DERIVATION OF THE TRANSPORT TERM
Let F be just convection term.
F 2t, x, u...., = - (v(*) x, u, v...)u
(23)
Because of the linearity of the inner product (-, 12, we can ignore other additive components of F in this derivation. In the rest of this section, we will write eg for o(k) (t, u, ...). We begin by plugging F into the message matrix from Eq. (14).
Mi -(V (uuu)),, (()))) n
(24)
Next, we split the domain into the mesh cells.
= --EAen (V (Duu)), ((0))
(25)
Now we apply the product rule.
=- AET. (ok Vuk + (V - 0k) uk. )((()))))))))))))))))))))))
(26)
We assume that the velocity field is iivergence-free, i.e. V V-ey = O.
ExET,
AET." (wk Vux, ( CH(A)
(27)
Now we expand u in the P1 basis (remember that we have encoded the data in the coefficients of u)
-'ser (Dk (uk-
N V j=1 yk
(0)) o (0)) CH(A)
(28)
and make use of the linearity of the inner product and gradient once more.
Σ SET.
Σ
N
Uk (vx Vipli), o ( (H(A)
(
(29)
j=1
Because the inner product is restricted to A, we can restrict the inner sum to vertices of A.
ESET, Σ L
yk ( (Ex Vi(li, )(()))
(30)
Finally, we factorize the inner products in the same way as we did it for the free-form term
ΣSET, Eames
Uk Un (Vili), (
(
(31)
and plug in the neural network go approximating the velocity field.
ESETT E
(1) (t,4) Uk Jo.A.A (Vi(0), (1) CH(A)
(32)
we stack this expression over all features and translate the result into the MPNN framework, we arrive at the transport term that we present as part of T-FEN in Section 3.1.
13

---

Published as a conference paper at ICLR 2022
C DATASETS
In this section, we give a brief overview over each of the datasets we used in Section 4 and how we pre-processed them.
Subsampling, Both real-world datasets that we use provide large amounts of dense data on 2D or 3D grid, which we subsample to simulate sparse measurements of a physical system and make the amount of data per timestep manageable. A non-trivial problem is the selection of sampling points and we decided to use the k-Medoids algorithm. This way we ensure that we use actual data points as sampling points and reach roughly uniform but still random cover of the domain. A uniform covering is especially important for the Black Sea dataset because of the non-convex shape of the domain and small bays the border that still want to cover.
C.I FLWW
(a) The magnitude the velocity field
The pressure field with the simulation mesh overlaid
Figure 7: A sample from the CylinderFlow dataset. The upper and lower boundaries are impenetrable walls. The left-hand boundary is an inflow with a constant velocity profile over time. The mesh resolution increases towards the boundaries at the top and bottom as well as in the immediate vicinity of the cylinder.
The CylinderFlow dataset is collection of simulated flows characterized by velocity and pressure fields released by Pfaff al. (2021). The domain two-dimensional channel with impenetrable boundaries at the top and bottom and circular obstacle of varying size and position that blocks the flow. On the left-hand boundary, the velocity field fixed to inflow profile that is constant over time but varies between samples. Depending on the obstacle, inflow profile, and flow parameters the resulting velocity and pressure fields exhibit distinct behaviors, varying between stable flows and oscillations. Pfaff et al. used the COMSOL solver to generate the data by solving the incompressible Navier-Stokes equations. See Fig. 7 for a visualization of the velocity and pressure fields as well as the simulation mesh.
The training set contains 1000 simulations, while validation and test set have 100 each. Every sequence consists of 600 steps of the flow sampled at time resolution of At 0.01s. Spatially, the data is discretized on meshes with about 1800 nodes. special property of this dataset is that each sequence has different mesh, because each mesh has notch representing the cylinder obstacle which varies size and position.
To enforce the boundary conditions, fix the velocity fields for all methods the boundaries except for the outflow, similar to Pfaff et al.. On the boundaries of the channel as well as the cylinder boundary, we fix the velocities zero. The velocity inflow profile stays constant over time and
14

---

Published as a conference paper at ICLR 2022
hold it fixed during training and inference. The pressure, however, is unconstrained throughout the domain and can also vary on the boundaries.
Due to the size of the dataset, we select only the first 100 sequences for training and train only subsequences starting at a multiple of 10 steps. in effect climinating overlapping training sequences and reducing the size of the training dataset by another factor of 10. For the evaluation, also choose the subsequences to test on in the same way, so every subsequence of length 10 starting at time step O, 10 and so on.
We normalize each feature by its mean and standard deviation across all training sequences.
C.2 BLACK SEA
The Black Sea data a reanalysis dataset of daily and monthly mean fields in the Black Sea of variables such as temperature, flow velocity, and salinity and 31 depth levels from 01/01/1992. The whole dataset is available as the Copernicus Marine Service under custom license" granting permission freely use the data in research under the condition of crediting the Copernicus Marine Environment Monitoring Service
1DXY__07_oa_ 20_00007002220000000 product from the EU
The spatial resolution of the raw data is 1/27 1/368 on a regular mesh of size 395x215 covering the Black Sea. Of the resulting 84,925 grid points, about 40,000 contain actual sea data and we subsample these valid points to create our sparse datasets. We mesh the chosen sample points with the Delaunay algorithm and filter acute mesh cells on the boundary as described in Appendix H. Additionally, we remove mesh cells that cover more land than water to avoid connecting nodes that are close by air distance but far apart through the sea.
We use daily data from 01/01/2012 to 31/12/2017 for training, 01/01/2018 to 31/12/2018 for validation and 01/01/2019 to 31/12/2019 for testing. As features we choose the zonal velocities and the water temperature. Both features are taken at a depth of 12.54m because the dynamics at the surface are strongly influenced by wind, sun, and cloud cover, which are not part of the data and thus make the dynamics non-deterministic from the perspective of the model.
On this dataset, we normalize the node positions as well as the features. The node positions are just normalized by their mean and average standard deviation in latitudinal and longitudinal direction over the training set as are the zonal velocity features. For the temperature, deviation grouped by calendar day over the training date range for the normalization, e.g. all January 1tss over the years 2012 to 2017 are grouped together, because the temperature exhibits a clear yearly cyclicity.
chose mean and standard
C.3 SCALARFLOW
ScalarFlow is a dataset of 104 smoke plumes published by Eckert et al. (2019) under the Creative Commons license CC-BY-NC-SA 4.0 (Eckert et al., 2019). Eckert et al. created a controlled environment in which fog from a fog machine rises in an air column over a heating element. They captured the resulting flow on video with calibrated multi-camera setup and reconstructed the velocity and density fields of each recording on a dense 100x178x (0xx178xx00 grid for 150 timesteps with their proposed method. At each grid point, the dataset provides 3-dimensional velocity vectors and the current fog density. The resulting dataset represents high-resolution measurements of dynamical physical system that is mostly driven by buoyancy and transport processes but also evokes turbulence. All plots of this dataset show the fog density.
Of the 104 recordings we use the first 64 for training, the next 20 for validation and the final 20 as a test set. Before subsampling the grid points, reduce the data to 2D by averaging over the depth (Kohl et al., 2020) and restrict the data to the central 60 points in x-direction and bottom 140 points in y-direction, because the fog does not reach points outside of that central box in almost any of the recordings. We normalize the coordinates of the sample points by their mean and standard deviation as well the features.
_00_0070004
15

---

Published as a conference paper at ICLR 2022
D MODEL CONFIGURATIONS
We trained all models on all datasets for at most 50 epochs or 24 hours. Most trainings have been stopped early after improvement in the validation score has been observed for 5 epochs the ODE-based models ran out of memory due to excessive memory requirements, when the learned dynamics required too many function evaluations to solve.
The three datasets present us with different circumstances regarding the basic properties of optimal dynamics that would model them well. Both Black Sea and ScalarFlow are best served by a non- stationary model, because it is reasonable to assume that the dynamics of the sea depend on the position and ScalarFlow has an inlet the bottom with fixed position. On top of that, we further expect that the BlackSea dataset has a time dependence exceeding the global, yearly fluctuations that we subtract with our normalization scheme. Therefore, we make all models non-stationary and non-autonomous Black Sea and just non-stationary on ScalarFlow. On CylinderFlow, the models are autonomous and stationary.
The mechanism to make FEN and T-FEN non-stationary and non-autonomous is described in Section 3.1. For the baselines, we concatenate the time stamps and the node positions respectively to the node features.
ODE-based models such as FEN and CT-MPNN require many function evaluations when their dynamics change abruptly. One possible trigger for that be jump in their inputs. This would, for example, occur if we would represent the time of the year on Black Sea as a number between O and to capture the yearly cyclicity of the dynamics. Then we would get discontinuity in the input at the turn of the year. To avoid this jump, we instead embed the time feature on Black Sea two-dimensionally (sin()),cos()) where i is a map from the current time within a year to []
For GWN, we use a batch size 6 on all datasets and a batch size of with PA-DGN. Due to an implementation detail of forchdiffes as of the submission of this paper, batched solving of ODEs with adaptive-step solvers can introduce interference between independent ODEs, because internal error estimates in the solver are computed across ODEs and the dynamics are stepped jointly. Therefore, we train and evaluate the ODE-based models with batch size of 1.
D.1 FEN & T-FEN
Both fo and go are MLPs with 4 hidden layers and tanh Iinaaiiaassies.. In FEN the layers of fo have a width of 128. For T-FEN we chose to reduce the width of the hidden layers of both fo and go to 96, so that both FEN and T-FEN have a comparable number of parameters and we can trace back any diffcrence in performance to difference model structure and not capacity. See Table for the parameter counts. We solve the feature dynamics with dopri5 solver with an absolute tolerance and relative tolerance 10-6 on CylinderFlow and Black Sea, and 10-7 on ScalarFlow, On initialization, set the last weights and bias of the last layers of both MLPs to 0 to start with the constant dynamics as a reasonable initialization for training. For parameter learning, we use the Adam optimizer with learning rate of 10-3 (Kingma & Ba, 2015).
D.2 BASELINES
CT-MPNN. Following lakovl al. (2021), our CT-
Table 2: Parameter counts for all models.
MPNN implementation models both the message func- tion and the update function as MLPs with tanh non- linearities. Contrary to the default configuration used by lakovi et (2021), we increased the widths of the MLPs from 60 to 128 and the message dimension from 40 to 128 to make the model capacity comparable to the other baselines. Instead of an implicit solver we used the same dopri5 solver as for FEN and T-FEN with an absolute tolerance of 10 and relative tolerance of 10 -7 This both increased the models performance and avoided frequent convergence issues in the step function of the -aahoohh solver. Finally, exchanged the RPROP opti-
Black Sea
ScalarFlow
PA-DGN
269140
269140
GWN
459750
466160
CT-MPNN
134403
134532
FEN
53129
53772
T-FEN
60975
61844
16

---

Published as a conference paper at ICLR 2022
mizer with standard Adam optimizer and we compute the gradients the same in FEN and T-FEN by backpropagation through the ODE solver.
We apply the same zero-initialization and increasing training sequence lengths, that use for our models, to this baseline.
PA-DGN. For PA-DGN choose two graph network layers of two layers each with hidden size of 64. We condition the model on 5 timesteps. To represent the structure of the domain, we pass the 4-nearest neighbor graph of the nodes to the model following authors of the method (Seo et al., 2020).
GWN. For GWN we use the improved configuration deseribed by Shleifer et al. (2019) with 40 convolution channels, extra skip connections and four of GWN's WaveNet inspired blocks with two layers each. The model is conditioned on 12 time steps data and predicts 10 steps ahead. For the long-term extrapolation Appendix G, we apply the model auto-regressively to get predictions for more than 10 steps. We do not learn an adaptíve adjacency matrix because that comes with quadratic memory cost in the number of nodes and Wu al. have shown that GWN only performs 1.2 % worse on METR-LA without the adaptive adjacency matrix. For the computation of the forward and backward transition matrices, we pass in the adjacency matrix derived from the mesh of the domain where the edges are weighted by the Gaussian kernel
[+(0)) - -0(i1112
e exp
(33)
where a is the standard deviation of ||x(00) x
E TRAINING
We train our models with the Adam optimizer (Kingma & Ba, 2015) by minimizing the multi-step prediction error in terms of the mean L1 distance between predicted and actual node features
1 NT
C(ú. y)
l (300)
(1:00 ll1
(34)
where (( is the prediction time ty. To accumulate gradients, backpropagate through the operations of the ODE solver, also known as discretize-ther seetiize-then-opee then-optimize as it improves training time significantly compared to backpropagation with the adjoint equation (Onken & Ruthotto, 2020).
We observed that training our models to predict complete sequences from the beginning can get them stuck in local minima that take many epochs to escape. To stabilize training, we begin by training on subsequences of length 3 and iteratively increase s by 1 per epoch until it reaches the maximum training sequence length.
F
For our super-resolution experiment, we have subsampled both datasets 3 times for each fixed number of nodes and evaluated each of the 10 trained models per model class against each of them. Therefore, the means (denoted by the markers) and standard deviations (denoted by the error bars) Fig. 3 and Fig. 8 have been computed over 30 evaluations.
Inspecting the full results in Fig. 8 reveals two things. While T-FEN performs best in super-resolution on both datasets, GWN and PA-DGN also only deteriorate slightly on finer meshes. However, this offset by the fact that these models did not achieve such a good to begin with and are still out-performed by T-FEN on every resolution. CT-MPNN achieved a lower prediction error than T-FEN on the original mesh granularity but generalizes markedly worse than both FEN and T-FEN to finer meshes. We suppose that this is due to the fact that CT-MPNN uses general MPNN internally, while FEN and T-FEN have an inductive bias towards producing physically sensible predictions because of their structural connection to FEM and PDEs.
17

---

Published as a conference paper at ICLR 2022
Black Sea
ScalarFlow
PADGN
MAB
GWN
CT-MPNN
FEN
T-FEN
0.09
2000
3000
4000
1000
3000
Nodes
Nodes
Figure 8: All models were trained on 1000 node subsamples of the data for 10-step prediction and then evaluated the same task but increasingly fine meshings of the domain.
Black Sea
ScalarFlow
1000 nodes
5000 nodes
Figure 9: A visual comparison of the data when subsampled at 1000 nodes and 5000 nodes.
G LONG TERM EXTRAPOLATION
In Fig. 10 we show 60-step extrapolations of the all models trained on 10-step forecasting. The first thing note is that PA-DGN struggles with the directionality inherent in the ScalarFlow dataset. While its prediction evolves over time, it exhibits a mostly diffusive dynamic.
From visual standpoint, GWN produces the most realistic looking forecasts with small-scale details even after many timesteps, especially compared to the ODE-based models that beat GWN in short-term forecasting. We
reasons for that. First, ODE-based models that learn the data dynamics directly, i.e. all three of FEN, T-FEN, and CT-MPNN, have to work with tight information bottleneck. Because these models carry information through time only the features themselves and have no latent state, they can condition their prediction only on a single time step and this bottleneck makes difficult to conserve fine details over time. GWN on the other hand predicts 10 timesteps at time and conditions each forecast the past 12 timesteps. In this way, GWN can extract more information from the input data and also carry that information forward 10 steps at time whereas ODE-based models have conserve the information through roughly 100 solver steps to reach f10. see the number of function evaluations in Table
Second, FEN, T-FEN, and CT-MPNN ODE models and such biased towards smooth predictions over time and over long time frames this leads to a smoothing in space. Yet, on short time frames, before smoothing sets these models are able model these physical processes more
18

---

Published as a conference paper at ICLR 2022
accurately, because ODEs are natural models for these data. For discussion of the differences between the FEN, T-FEN, and CT-MPPN predictions, see Section 4.
Data
FEN
T-FEN
CT-MPPNN
GWN
PADGN
to + 805ms
to + 536ms
to + 268ms
to + 13ms
Figure 10: Predictions of all models for rising smoke plume after 20, 40 and 60 time steps.
H FILTERING HIGHLY ACUTE CELLS ON THE BOUNDARY
Delaunay triangulation of set of points can produce very long cells on the boundary as in Fig. 11 when several points close to the boundary are almost co-linear, This can just occur in the data but also happen during generation of synthetic data. If some boundary nodes in synthetic data are exactly co-linear such as in grid domains, they might lose the exact co-linearity due to the inexactness of floating point numbers and operations if we, for example, rotate the domain. To ensure that these artifacts do not impact our models, we pre-process the triangulation with Algorithm and an angle threshold of 10x
180
19

---

Published as a conference paper at ICLR 2022
Algorithm 1 Filtering cells with small interior boundary-adjacent angles from a triangulation given threshold angle and
dimensional boundary faces
1: function FILTER FILTER Stainnaaaannnnnnnnaaaaaanon T, threshold e > 0, boundary faces B) for all A
2:
T do
> Count number of faces the boundary for each cell
3:
Â c A, | = E B)]
4:
F+ { TTcs = =}
5:
while [F] do
6:
A < pop(FF
7:
if # then
Check that A is not edge corner of the domain
8:
continue
9:
(Á.E) < (A)it(A)
Split boundary face and interior node
BONNDOUUDDAYYANGLLL,,, { < T + [4]
then
11:
Remove cell Add newly become boundary faces to B
12:
for all Á C Á ( d
B BU (A U (E]}
for all A' T A A'
15:
+ ca +1
16:
FU [N}
Add cells that are now boundary cells to J
17:
return T
18:
function MINHOOUDDARYANGLERS....
19:
B + project(5, Y イイズ π
Project & onto the hyperplane given by Á
for all A E Â do
22:
+ min(1 y,, AC)
Find the angle between the line segments AB and Ag
23:
return 7
20

---

Published as a conference paper at ICLR 2022
(a) The inside angles of the shaded boundary cell below the threshold and is removed.
(b) The removal has made this (previously) interior cell into boundary cell that itself falls below the angle threshold and gets removed as well.
Figure 11: Delaunay triangulation can produce elongated, splinter-like cells when interior nodes are close to the boundary of the convex hull as in this example. We post-process the Delaunay triangulations with Algorithm filter these cells out.
21


Reviews:
[{"review_id": "Nd6R_4_o-Gq", "reviewer": "Reviewer_4K4v", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "This manuscript proposes a new graph neural net (GNN) method to learn the dynamics of a spatio-temporal PDE-driven dynamical system directly from data. The authors propose to do that using the finite element method (FEM). The proposed method builds on using: basis function approximation for the (unknown) field u, Galerkin method with the assumption that discrepancy between the dynamics F and basis function approximation is orthogonal to (finite) basis functions, method of lines, and message passing GNN as a proxy for the dynamics. The use of linear interpolation allows to express the time derivative of $Y$ as a solution of a system of linear equations, which is further approximated to gain computational additional efficiency. Authors also propose a method to incorporate inductive bias into model learning for models that are assumed to contain a convection component. Overall the proposed method is well-motivated, and for the most part the description is clear. To my knowledge, the proposed method is novel and contains some methodologically new ideas, and the performance seems to be on par with previous methods that learn free-form dynamics, and shows an improvement for models that contain a convection component when such prior knowledge is utilised in the model training. Authors could address and/or clarify the following aspects:\n1. I understand that a piece-wise linear basis simplifies computational complexity by making some of the computational steps straightforward, but on the other hand the selected basis is the simplest and obviously not optimal from approximation accuracy point of view. Can the method be extended to other basis? For example, if we knew that the dynamics $F$ contains a diffusion term $\\nabla^2u$ we would not be able to introduce it since the second derivative of PWL functions is zero everywhere. I think the discussion of possible limitations of the PWL basis and of possible extensions to higher-order bases is missing.\n2. It seems that the measurement from the initial time point is used as the initial state? Why not introduce a separate free parameter for the initial state? \n3. As described towards the end of the manuscript, the system is initialed with an initial state and then the PDE dynamics define how the system evolves over time. However, below eq. (12) it is noted that \"where $u$ encodes the data $Y$ at time $t$ in its coefficients.\" Is $u$ defined based on the dynamically evolving system state or by data? \n4. Paragraph \"Network architecture\" on page: scimitar-learn is used to compute inner product between basis functions. Provide a brief description of how the computation is done. \n5. End of page 6: gradients are back propagated through an ODE solver. Why not use adjoint method (possibly with checkpoints) as proposed in previous work and implemented e.g. in comparison methods? Discrete back-propagation may not scale to longer sequences. Is this the reason why data trajectories are shortened? \n6. Table 1: split the methods into two groups. Group 1 should includes PA-DGN, GWN, CT-MPNN and FEN, which do not assume prior knowledge about system dynamics. FEN performs similarly with CT-MPNN on ScalarFlow and perhaps slightly better on Black Sea data (MAE values within one std). Group 2 should include only T-FEN that is specifically designed to learn convection systems, and thus provides a small performance improvement. \n7) In experiments, the authors make their models time- and position-dependent while the strongest baseline models (CT-MPNN) does not utilize neither time nor positions. That makes it hard to tell whether the improvements in performance of FEN and T-FEN are due to the models's structure and inductive biases or due to time and position dependencies. Authors should provide an ablation study to address this.\n8. Provide additional comparisons on systems with larger variety of dynamics using simulated data (for which the ground truth is known) to better understand when FEN performs better/worse than comparison methods. \n9. Datasets are subsampled to 1000 spatial points. Provide results for smaller and also larger spatial grids to demonstrate the argument that \"approximation becomes arbitrarily good as the mesh resolution increases\".  ", "summary_of_the_paper": "This manuscript proposes a new graph neural net (GNN) method to learn the dynamics of a spatio-temporal PDE-driven dynamical system directly from data.", "summary_of_the_review": "To my knowledge, the proposed method is novel and contains some methodologically new ideas, and the performance seems to be on par with previous methods that learn free-form dynamics, and shows an improvement for models that contain a convection component when such prior knowledge is utilised in the model training. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "NNiEtvWHcAh", "reviewer": "Reviewer_Zk77", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "# Strenghts\n\n- The paper looks at a very realistic setting for learning PDEs from data: a finite number of samples and (partially) unknown true dynamics. Tackling these problems concomitantly is of great practical importance and this makes the proposed method relevant for practical applications. \n- The introduction does a good job motivating the work and pinpointing the main challenges of learning dynamics from data. \n- As someone with little experience with the finite element method, Section 2.1 does a great job explaining the required background in the right amount of detail for understanding the paper.\n- The connection that the paper makes between the finite element method and message passing neural networks is interesting and, to the best of my knowledge, original. \n- The authors show how inductive biases can be added to the model by using a certain prior over the structure of the function $F$. \n- I like that the paper focuses on real-world datasets. Also, the datasets themselves are extremely interesting to visualize and make the paper more interesting. \n- Figure 3 provides a very insightful qualitative understanding of the proposed model compared to the baseline. \n- I am glad that the paper includes a super-resolution experiment. Often, models that work with a discretised space can be very sensitive to changes in the resolution of the mesh. Figure 4 shows that the proposed model is relatively robust to changes in the number of triangles. \n- The authors show (in a relatively specific setting) that factorizing the dynamics achieves a disentanglement effect which allows some degree of interpretability of the model. \n\n# Weaknesses \n\n- While I appreciate the focus on real-world datasets, a synthetic experiment where the model could have been evaluated in a more systematic way would have been useful. \n- The trick to stabilize training described at the end of Section 3 is slightly peculiar. How important is this trick? Do the authors have any results for when this trick is not used? Could this trick improve the performance of the baselines as well? \n- The importance of the approximation from Equation (13) is not studied. Perhaps that is something that could have been tried in the more controllable synthetic setting I was suggesting above. In general, I would be interested to know what are the costs of this approximation and if more advanced approximations might be worth being considered to boost performance.\n- Minor suggestion: The paper frames the model as a hypergraph neural network. However, the authors might want to be aware that there is a recent line of work developing simplicial and cell complex neural networks: https://arxiv.org/abs/2103.03212 (ICML 2021), https://arxiv.org/abs/2106.12575 (NeurIPS 2021), https://arxiv.org/abs/2010.03633. Since the model learns a function over the 2-simplices in the simplicial complex, the model is probably more accurately described as a type of simplicial neural network.   ", "summary_of_the_paper": "The paper proposes a graph / simplicial neural network based on the finite element method for learning dynamics from data when only a finite number of samples exist and the true dynamics are not known or only partially known.  ", "summary_of_the_review": "The weaknesses reported above are relatively minor and far outweighed by the strengths of the paper. Therefore, I recommend the paper for acceptance. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "HdrzW1ut77", "reviewer": "Reviewer_EPhP", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Strengths:\n\n- The paper is well written and well presented. On the whole it was relatively easy to understand and the diagrams definitely contribute to that.\n- The proposed model is well motivated and nicely extends ideas such as Graph Neural ODEs to the PDE domain.\n- The paper considers the issues with a naive implementation in depth, (which is that the model would be slow). It then provides the solution to this.\n- Rigorous tests on the proposed model are carried out against robust baselines.\n- There is an extensive review of related work.\n- A lot of effort has been put in to make the results reproducible, including all experimental details and code/datasets to come.\n\nWeaknesses/Questions:\n\n- My main concern is that training has been carried out over 10 time steps, was this a hyperparameter that was tuned? I agree that the correlation over say 30 steps will be minimal, however my understanding is that all models learn an update based on the current state and potentially a few steps in the past. Would it not still be possible to train over more time steps? Could T-FEN extrapolate better than GWN if it is trained on a larger time range?\n- One of the proposed reasons GWN outperforms T-FEN at extrapolation is that it can use the past time-steps as input. This could be extended to T-FEN in the form of delay differential equations (https://arxiv.org/abs/2102.10801), where the dynamics takes the state at $(t)$ and the state at $(t-\\tau)$ as input. This would make the model use the past as well. Would it be possible to carry out an ablation study using this?\n- The paper could benefit from a discussion section. Saying what the model is good at and bad at. For example, we see that it performs very well in the times used to train, but is not as good at extrapolation.\n- In table 1, for the ODE based method, the number of function evaluations are provided to show T-FEN is faster than FEN and CT-MPNN. Is it possible to time the evaluation to support this claim further? That way one could also compare to the other two baselines GWN and PA-DGN.\n- Is it possible to move the related work section, either to the introduction or to just before the conclusion? It breaks the flow.\n- What is the reason for using the $L_1$ loss over the $L_2$ loss? Could we expect better/worse results with $L_2$?\n- Could the model be improved even further by using more physically informed terms? For example, like those in electro-magnetism $\\dot{v}=qv \\times B$, where $B$ is learnt (https://arxiv.org/abs/2109.07359). Or, because the experiments use data of fluids, could it help to include Laplacian or curl terms in the dynamics, which appear in Navier-Stokes? Is it possible to carry out an ablation? \n- The mass matrix $A$ is approximately inverted by lumping the matrix, are there situations where this approximation could lead to errors?\n- Is it possible to extend this method where we expect higher order PDEs, that include terms such as $\\frac{\\partial^{2}u}{\\partial t\\partial x}$?\n- The model disentangles the dynamics into a convection term and the remainder. How does this relate to disentangling the dynamics of an Augmented Neural ODE (https://arxiv.org/abs/1904.01681) into a velocity and an acceleration (https://arxiv.org/abs/2006.07220), which can also aid the interpretability of the ODE?\n- I don’t entirely understand the Black Sea dataset. Is the mean temperature the mean temperature of the entire sea or of a region, say $1m^2$? If it is of the whole sea, does this not remove the spatial element of the task? Additionally are there any interesting effects appearing over long time periods due to global warming/concept drift? Given that the training regime is taken over 2012-2017 and the testing regime is in 2019?\n\nMinor Points:\n\n- The paragraph just above equation 15 has a small mistake. The manuscript says “This would prohibit training even with the adjoint equation (Chen et al., 2018)”, while the paragraph is about the speed of training. The benefit of the adjoint method is that it is memory efficient but slow compared to directly backpropagating through an ODE solver, which is fast but uses a lot of memory.\n- There is a typo at the bottom of page 8, “disappears Extrapolationat the”\n- There is a typo at the bottom of page 16, “ODE-base models” should be “ODE-based models”\n", "summary_of_the_paper": "This paper proposes a new model for learning partial differential equations from data. The PDE is first discretized then solved as an ODE. The dynamics function is learned with Message-Passing Neural Networks, where the function is split into a sum of physically informed terms. This splitting both improves model performance and makes the model more interpretable by disentangling the dynamics. The model is tested rigorously against multiple baseline models, and the results show the new model performs well.", "summary_of_the_review": "The paper is well written, the model (to my best knowledge) is novel, with the method building on existing work. The experiments rigorously test the model against the necessary baselines and information is given in the appendices on reproducing the results. Therefore I recommend acceptance, with a few clarifications to be made.\n\n**EDIT** I have increased my confidence score from 3 to 4 after my initial questions have been answered.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "4HoSyF0mYsy", "reviewer": "Reviewer_CGtm", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "The paper uses Message Passing Neural Networks to implement a Finite Element Method with learnable dynamics. It combines different models and techniques from the literature, but it clearly does that in a nontrivial way including modifications, therefore it is more than just derivative work.\n\nI find it valuable that the author also considers performance of implementation, and discusses consideration on GPU architecture related performance issues. Also the effort put into making the work reproducible adds value to the paper.\n\nThe method can clearly have practical value, and the discussion of the method is clear and quite detailed. I am not as familiar with the PDE/NN literature as for example the Neural ODE literature, therefore I cannot rule out entirely that something similar already exists. \n\nQuestions / actionable comments:\n\nI)\n\"The results show that FEN and T-FEN provide the smallest prediction error on both datasets with a further boost due to the separate transport term in T-FEN.\"\nI feel this is a bit too strong statement. While I tend to accept the paper with the results in Table 1 given the other clear benefits of the method (like robustness and interpretability) I am not comfortable with this statement. The table is based on 3 repeats, and have results for example:\n\nBlack Sea dataset  CT-MPNN 0.944 ± 0.003  vs.  FEN 0.938 ± 0.005, \n\nthese confidence intervals are clearly overlapping and using 3 samples I am a bit skeptical. This is one of the strongest gaps in the table.\nIn the ScalarFlow FEN is really on par with CT-MPNN. Did you used a statistical test to bold the results?\n\nAgain, I have no problem supporting the paper even if the method is on par with the best competitor, but the statement in the paper should be supported by the statistics. \n\nSimilarly, in the case of NFE CT-MPNN have approx +-60 deviation on the ScalarFlow dataset, I am not too comfortable to compare these numbers.\n\nII)\nWhat is the motivation of using L1 loss? instead of like MSE?\n\nIII) In Experiments/Multi-step Forecasting section:\nTime horizon choice is well motivated, but still feels arbitrary in some sense. Does the author see some way to formalize what we accept as \"meaningful dynamics\"?\nHow one should choose a comparison horizon for example. Can it be done at least approximately without domain knowledge on the system? \n\nIV) In Model/Network architecture section: \nOrder invariance of cell vertices is assured by ordering the nodes canonically. Would the author expect improvement if a permutation invariant network, like a Set Transformer be used here? Why or why not?\nIn some sense the message aggregation step, summation being permutation invariant, is a set network, but not a very expressive one. This  however makes the cell order invariant for a given node and not the other way around. \n", "summary_of_the_paper": "The author proposes a method for forecasting in Partial Differential Equations by coupling Finite Element Method on an arbitrary grid with the learning of the dynamics from data. For this purpose a variant of message passing based graph networks is used. It is show in the paper that it's possible to incorporate priors on the structure of the PDE that results in an interpretable solution. The model also show more stability to changes of the mesh structure in test time (like superresolution) and to extrapolation than competitors.\n", "summary_of_the_review": "The paper give valuable contribution, the method expected to be practical, robust, and in some cases interpretable. I find the statement on raw prediction error overly strong.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper introduces a graph neural network (GNN) based on the finite element method (FEM) for learning partial differential equations from data. The proposed finite element network is based on a piecewise linear function approximation and a message passing GNN for dynamics' prediction. The authors also propose a method to incorporate inductive bias when learning the dynamical model, e.g. including a convection component.

The paper received three clear accept and one weak accept recommendations. The reviewers discussed the possible extensions of the method, and also raise several concerns regarding experiments, e.g. the added value of a synthetic dataset, implementation tricks or hyper-parameter settings. The rebuttal did a good job in answering reviewers' concerns: after rebuttal, there was a consensus among reviewers to accept the paper.

The AC's own readings confirmed the reviewers' recommendations. The paper is well written and introduces solid contribution at the frontier of GNNs and finite elements methods, especially a pioneer graph-based model for spatio-temporal forecasting derived from FEM. Therefore, the AC recommends acceptance.
Decision: Accept (Spotlight)