Published as a conference paper at ICLR 2022

T OWARDS FAST S IMULATION OF
E NVIRONMENTAL F LUID M ECHANICS WITHM ULTI -S CALE G RAPH N EURAL N ETWORKS
Stathi Fotiadis
Department of Bioengineering
Imperial College London

Mario Lino
Department of Aeronautics
Imperial College London
mal1218@ic.ac.uk

Chris Cantwell
Department of Aeronautics
Imperial College London

Anil A. Bharath
Department of Bioengineering
Imperial College London

A BSTRACT
Numerical simulators are essential tools in the study of natural fluid-systems, but
their performance often limits application in practice. Recent machine-learning
approaches have demonstrated their ability to accelerate spatio-temporal predictions, although, with only moderate accuracy in comparison. Here we introduce
MultiScaleGNN, a novel multi-scale graph neural network model for learning toinfer unsteady continuum mechanics in problems encompassing a range of lengthscales and complex boundary geometries. We demonstrate this method on advection problems and incompressible fluid dynamics, both fundamental phenomena
in oceanic and atmospheric processes. Our results show good extrapolation to new
domain geometries and parameters for long-term temporal simulations. Simulations obtained with MultiScaleGNN are between two and four orders of magnitude
faster than those on which it was trained.

1

I NTRODUCTION

Forecasting the spatio-temporal mechanics of continuous systems is a common requirement in many
areas of science and engineering, including environmental fluid dynamics (Rubin, 2001). Physicalmodels for the transport and dispersion processes in natural fluid flows often consist on one or morepartial differential equations (PDEs), whose complexity may preclude their analytic solution (Kim
& Boysan, 1999). Numerical methods are well-established for approximating the solution of PDEs
with high accuracy, but they are computationally expensive (Karniadakis & Sherwin, 2013). Deeplearning techniques have been shown to be capable of accelerating physical simulations (Guo et al.,2016). Most of the recent work on deep learning to infer continuum physics has focused on theuse of convolutional neural networks (CNNs). In part, the success of CNNs for these problems
lies in their translation invariance and locality Goodfellow et al. (2016), which represent strong anddesirable inductive biases for continuum-mechanics models. However, CNNs constrain input andoutput fields to be defined on rectangular domains represented by regular grids, which is not suitable
for more complex domains. As with traditional numerical techniques, it is desirable to be able tovary the resolution in space, devoting more effort where the physics are challenging to resolve,
and less effort elsewhere. An alternative approach to applying deep learning to geometrically andtopologically complex domains is provided by graph neural networks (GNNs), which can also be
designed to satisfy spatial invariance and locality (Battaglia et al., 2018; Wu et al., 2020).
In this paper, we describe a novel approach to applying GNNs for accurately forecasting the evolution of physical systems in complex and irregular domains. We propose MultiScaleGNN, a multiscale GNN model, to forecast the spatio-temporal evolution of continuous systems discretised as
unstructured sets of nodes. Each scale processes the information at different resolutions, enabling
the network to more accurately and efficiently capture complex physical systems. We apply Mul1

Published as a conference paper at ICLR 2022

tiScaleGNN to simulate advection and incompressible fluid dynamics, fundamental processes inoceanic and atmospheric research. Importantly, MultiScaleGNN is independent of the spatial discretisation and the mean absolute error (MAE) decreases linearly as the distance between nodes isreduced, which allows for the application of adaptive re-meshing (see Figure 1b) MultiScaleGNN
simulations are between two and four orders of magnitude faster than the numerical solver used forgenerating the ground truth datasets, becoming a potential surrogate model for fast predictions.

(a) Advection(b) Incompressible flow

Figure 1: MultiScaleGNN forecasts continuum dynamics. We used it to simulate (a) advection [video] and,(b) the incompressible flow around a circular cylinder [video].

2

R ELATED WORK

During the last five years, most deep neural networks (DNNs) architectures used for predicting fluid
dynamics have included convolutional layers (Guo et al., 2016; Tompson et al., 2017; Lee & You,
2019; Kim et al., 2019; Wiewel et al., 2019; Fotiadis et al., 2020). These CNN-based solvers are
between one and four orders of magnitude faster than numerical solvers (Guo et al., 2016; Lino
et al., 2020), and some of them have shown good extrapolation to unseen domain geometries andinitial conditions (Thuerey et al., 2018; Lino et al., 2020). Recently, GNNs have been used to simulate the motion of discrete systems of solid particles (Battaglia et al., 2016; Chang et al., 2017) anddeformable solids and fluids discretised into Lagrangian (or free) particles (Li et al., 2019a; Mrowca
et al., 2018b; Sanchez-Gonzalez et al., 2020). Further research in this area introduced more general message-passing (MP) layers (Sanchez-Gonzalez et al., 2018; Li et al., 2019b; Mrowca et al.,2018a), high-order time-integration (Sanchez-Gonzalez et al., 2019) and hierarchical models (Liet al., 2019a; Mrowca et al., 2018b). To the best of our knowledge Alet, et al. (2019) Alet et al.
(2019) were the first to explore the use of GNNs to infer Eulerian mechanics by solving the Poisson PDE. However, their domains remained simple, used coarse spatial discretisations and did not
explore the generalisation of their model. More closely related to our work, Pfaff et al. (2021) proposed a mesh-based GNN to simulate continuum mechanics, although they did not consider the useof MP at multiple scales of resolution. Li et al. (2020b) and Liu et al. (2021) used multi-resolution
GNNs to infer steady solutions, but their pooling remained simple and did not explore extrapolation
to unseen domain geometries or physical parameters. Li et al. (2020a) also considered multi-scaleneural PDE modelling, with the mayor drawback that the spatial discretisations are constrained touniform grids.

3

M ODEL

3.1

M ODEL DEFINITION2
For a PDE ∂u∂t = F(u) on a spatial domain D ⊂ R , MultiScaleGNN infers the temporal evolution
1of u(t, x) at a finite set of nodes V , with coordinates x1i ∈ D. Given an input field u(t0 , xV 1 ), attime t = t0 and at the V 1 nodes, a single evaluation of MultiScaleGNN returns u(t0 + dt, xV 1 ),where dt is a fixed time-step size. Each time-step is performed by applying MP layers in L graphs
and between them, as illustrated in Figure 2. The high-resolution graph G1 consists of the set ofnodes V 1 and a set of directed edges E 1 connecting these nodes. In a complete graph, there exist

2

Published as a conference paper at ICLR 2022

Figure 2: MultiScaleGNN architecture. G1 is the high-resolution graph, G2 and GL are low-resolution graphs.
|V 1 |(|V 1 |−1) edges, this would result in a high computational cost to MP. Instead, MultiScaleGNN
connects each node in V 1 (receiver node) to its k-closest nodes (k sender nodes) using a k-nearest
neighbours (k-NN) algorithm. This guarantees that, in MP layers, each node is receiving information from exactly k nodes, in contrast to Pfaff et al. (2021) and Sanchez-Gonzalez et al. (2018).
This uniformity can ease the learning and improve the accuracy of the predictions. The attributes
assigned to each node vi1 are the concatenation of u(t0 , xi ), pi and Ωi , where pi is a set of physicalparameters at xi (such as the Reynolds number in fluid dynamics) and Ωi = 1 for nodes locatedon Dirichlet boundaries and zero elsewhere. Each edge attribute e1k is assigned the relative position
between sender node sk and receiver node rk .
Node attributes and edge attributes are encoded through two independent multi-layer perceptrons(MLPs). A MP layer applied to G1 propagates the nodal and edge information only locally between
adjacent nodes. Nevertheless, most continuous physical-systems require this propagation at larger
scales, or even globally. To handle this, MultiScaleGNN processes the information at L scales,
creating a graph for each level and propagating information between them in each pass. The lowerresolution graphs (G2 , G3 , . . . , GL ; with |V1 | > |V2 | > · · · > |VL |) possess fewer nodes and edges,
and hence, a single MP layer can propagate features over longer distances more efficiently. As
depicted in Figure 2, the information is first diffused and processed in the high-resolution graph
G1 through M1 MP layers. It is then passed to G2 through a downward MP (DownMP) layer. In
G2 the attributes are again processed through M2 MP layers and a DownMP layer to G3 . This
process is repeated L − 1 times. The lowest resolution attributes (stored in GL ) are then passed
back to the scale immediately above through an upward message passing (UpMP) layer. Attributes
are again successively passed through Ml MP layers at scale l and an UpMP layer from scale l toscale l − 1 until the information is ultimately processed in G1 . Finally, a MLP decodes the nodal
information to return the predicted field at time t0 +dt at the V 1 nodes. To apply MP in the L graphs,
MultiScaleGNN uses the MP layer introduced by Sanchez-Gonzalez et al. (2018) and Battaglia et al.
(2018), with the mean as the aggregation function.
3.2

M ULTI - SCALE GRAPHS

Each graph Gl = (V l , E l ) with l = 2, 3, . . . , L; is obtained from graph Gl−1 by first dividing D
into a regular grid with cell size dlx × dly (Figure 3a exemplifies this for l = 2). For each cell,provided that there is at least one node from V l−1 in it, a node is added to V l . The nodes from V l−1
in a given cell i and the node from V l on the same cell are denoted as child nodes, Ch(i) ⊂ V l−1 ,and parent node, i ∈ V l , respectively. The coordinates of each parent node is the mean positionof its children nodes. Each edge k ∈ E l connects sender node sk ∈ V l to receiver node rk ∈ V l ,provided that there exists at least one edge from Ch(sk ) to Ch(rk ) in E l−1 . Edge k is assigned themean edge attribute of the set of edges going from Ch(sk ) to Ch(rk ).
Downward message-passing (DownMP). To perform MP from Gl−1 to Gl (see Figure 3b), a setof directed edges, E l−1,l , is created. Each edge k ∈ E l−1,l connects node sk ∈ V l−1 to its parent
node rk ∈ V l , with edge attributes assigned as the relative position between child and parent nodes.
A DownMP layer applies a common edge-update function, f l−1,l , to edge k and node sk . It thenassigns to each node attribute vil the mean updated attribute of all the edges in E l−1,l arriving tonode i, i.e.,
3

Published as a conference paper at ICLR 2022(a) Building V 2 from V 1(b) DownMP and UpMP layers

Figure 3: (a) V 2 is obtained from V 1 by partitioning D using a grid with cell size d2x × d2y . (b) DownMP and
UpMP diagrams.

vil =

X
1
f l−1,l ([el−1,l, vsl−1]), ∀i ∈ 1, . . . , |V l |.
kk
|Ch(i)|(1)

k:rk =i

Upward message-passing (UpMP). To pass and process the node attributes from Gl+1 to Gl ,MultiScaleGNN defines a set of directed edges, E l+1,l . These edges are the same as in E l,l+1 , butwith opposite direction. An UpMP layer applies a common edge-update function, f l+1,l , to eachedge k ∈ E l+1,l and both its sender (at scale l + 1) and receiver (at scale l) nodes, directly updating
the node attributes in Gl , i.e.,
vrl k = f l+1,l ([el+1,l, vsl+1, vrl k ]), ∀k ∈ E l+1,l .
k
k

(2)

UpMP layers leave the edge attributes of Gl unaltered. To model functions f l−1,l and f l+1,l , weuse MLPs.

4

T RAINING DATASETS

Datasets AdvBox and AdvInBox, both used simultaneously for training, contain simulations ofa scalar field advected under a uniform velocity field on a square domain ([0, 1] × [0, 1]) and arectangular domain ([0, 1] × [0, 0.5]) respectively. AdvBox domains have periodic conditions on
all four boundaries, whereas AdvInBox domains have upper and lower periodic boundaries, aprescribed Dirichlet condition on the left boundary, and a zero-Neumann condition on the rightboundary. The initial states at t0 are derived from two-dimensional truncated Fourier series withrandom coefficients and a random number of terms. For advection models, u(t, xi ) ∈ R is theadvected field and pi ∈ R2 are the two components of the velocity field at xi .
Dataset NS contains simulations of the periodic vortex shedding around a circular cylinder at
Reynolds numbers between 500 and 1000. The upper and lower boundaries are periodic, and the vertical distance between cylinders is randomly sampled between 4 and 6. Each domain is discretised
into approximately 7000 nodes. For flow models, u(t, xi ) ∈ R3 contains the velocity and pressure
fields and pi ∈ R is the Reynolds number. Further details of the training and testing datasets areincluded in Appendix A.

5

R ESULTS AND DISCUSSION

We first consider the MultiScaleGNN model trained to infer advection. Despite MultiScaleGNN
being trained on square and rectangular domains and uniform velocity fields, it generalises to complex domains and non-uniform velocity fields. As an example of a closed domain, we consider the
Taylor-Couette flow in Figure 4, where after 49 time-steps a MultiScaleGNN with L = 1 maintains
high accuracy in transporting both the lower and the higher frequencies generated due to the shear
flow. We also evaluated the predictions of MultiScaleGNN on open domains containing obstaclesof different shapes (circles, squares, ellipses and closed spline curves). Figure 1 shows the predictions obtained with a MultiScaleGNN with L = 2 for a field advected around a circular cylinderand upper-lower and left-right periodic boundaries. MultiScaleGNN was also trained to simulate
4

Published as a conference paper at ICLR 2022

unsteady incompressible fluid dynamics in a range of Reynolds numbers between 500 and 1000. It
shows very good interpolation to unseen Reynolds numbers within this range. For instance, Figure
1b shows the horizontal velocity field predicted by MultiScaleGNN with L = 4 after 99 time-steps
for Re = 700. The MAE increases for Reynolds numbers lower than 400 and higher than 1200, but
the predictions remain visually realistic in the range 300 to 1500, indicating that the network haslearnt some fundamental aspects of the fluid dynamics. 1
We evaluate the accuracy of MultiScaleGNN with L = 1, 2, 3
and 4; the architectural details of each model are includedin Appendix B. Tables 1 and 2 collect the MAE for the lasttime-point and the mean of all the time-points on the testingdatasets. Incompressible fluids have a global behaviour, andthe addition of coarser layers helps the network to learn thischaracteristic and achieve significantly lower errors. Hence,for the N-S testing dataset there is a clear benefit from increasing the number of scales. In contrast, for the advection
datasets, the lowest MAE values are obtained for L = 1 or 2,
since in advection the information is propagated only locally.
As a comparison to Pfaff et al. (2021), a GNN with 16 sequential MP-layers (GN-Blocks) results in a MAE of 5.852 × 10−2
on our NSMidRe dataset; whereas MultiScaleGNN with the Figure 4: Advection of a scalar fieldsame number and type of MP layers, but distributed among 3 in a Taylor-Couette flow [video].
scales, results in a lower MAE of 3.081 × 10−2 . A comparisonof our coarsening/pooling algorithm to Liu et al. (2021) is included in Appendix C.
Table 1: MAE ×10−2 on the advection testing datasets for MultiScaleGNN models with L = 1, 2, 3, 4
Datasets

L=1
Step 49
All

L=2
Step 49
All

L=3
Step 49
All

L=4
Step 49
All

AdvTaylor
AdvCircle
AdvCircleAng
AdvSquare
AdvEllipseH
AdvEllipseV
AdvSplines
AdvInCir

6.940
4.057
3.870
4.250
4.462
4.265
4.484
11.7

7.914
3.690
3.962
4.175
4.328
4.132
4.426
25.369

8.037
3.911
4.300
4.420
4.567
4.297
4.636
27.981

8.790
3.827
4.428
4.279
4.449
4.207
4.58
27.379

3.195
2.112
2.030
2.178
2.241
2.236
2.293
7.101

3.676
1.817
1.890
1.991
2.017
2.014
2.088
18.548

3.739
1.916
2.074
2.107
2.136
2.100
2.199
22.626

4.050
1.995
2.249
2.141
2.179
2.22
2.300
22.034

Table 2: MAE ×10−2 on the N-S testing datasets for MultiScaleGNN models with L = 1, 2, 3, 4
Datasets

L=1
Step 99
All

L=2
Step 99
All

L=3
Step 99
All

L=4
Step 99
All

NSMidReNSLowRe
NSHighRe

9.765
9.43
9.487

4.759
12.3467.879

3.851
11.707
6.98

3.456
7.338
5.826

6

6.108
7.327
10.598

3.663
8.211
9.002

3.081
8.203
7.096

2.825
5.532
5.871

C ONCLUSION

MultiScaleGNN is a novel multi-scale GNN model for inferring mechanics on continuous systems
discretised into unstructured sets of nodes. Unstructured discretisations allow complex domains tobe accurately represented and the node count to be adjusted over space. Multiple coarser levels allow
high and low-resolution mechanics to be efficiently captured. In global and local problems, such as
incompressible fluid dynamics, the coarser graphs are particularly advantageous, since they enable1

Animations with the ground truth and best-model predictions can be found here.

5

Published as a conference paper at ICLR 2022

global characteristics to be learnt. MultiScaleGNN interpolates to unseen spatial discretisations of
the physical domains, allowing it to adopt efficient discretisations and to dynamically and locallymodify them to further improve the accuracy. MultiScaleGNN also generalises to advection oncomplex domains and velocity fields and it interpolates and extrapolates to unseen Reynolds numbers in fluid dynamics. Inference is between two and four orders of magnitude faster than with thehigh-order solver used for generating the training datasets. This work is a significant advancementin the design of flexible, accurate and efficient neural simulators for fluid dynamics and in generalfor continuum mechanics.

R EFERENCES
Ferran Alet, Adarsh Keshav Jeewajee, Maria Bauza Villalonga, Alberto Rodriguez, Tomas LozanoPerez, and Leslie Kaelbling. Graph element networks: adaptive, structured computation andmemory. In International Conference on Machine Learning, pp. 212–222. PMLR, 2019.
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint
arXiv:1607.06450, 2016.
Peter W Battaglia, Razvan Pascanu, Matthew Lai, Danilo Rezende, and Koray Kavukcuoglu. Interaction networks for learning about objects, relations and physics. In 38th Neural Information
Processing Systems, NeurIPS 2016, 2016.
Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi,
Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al.
Relational inductive biases, deep learning, and graph networks. arXiv preprint arXiv:1806.01261,
2018.
Filipe de Avila Belbute-Peres, Thomas Economon, and Zico Kolter. Combining differentiable pde
solvers and graph neural networks for fluid flow prediction. In International Conference on Machine Learning, pp. 2402–2411. PMLR, 2020.
Chris D Cantwell, David Moxey, Andrew Comerford, Alessandro Bolis, Gabriele Rocco, Gianmarco
Mengaldo, Daniele De Grazia, Sergey Yakovlev, J-E Lombard, Dirk Ekelschot, et al. Nektar++:
An open-source spectral/hp element framework. Computer physics communications, 192:205–
219, 2015.
Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum. A compositional
object-based approach to learning physical dynamics. In 5th International Conference on Learning Representations, ICLR 2017, 2017.
Stathi Fotiadis, Eduardo Pignatelli, Mario Lino, Chris Cantwell, Amos Storkey, and Anil A. Bharath.
Comparing recurrent and convolutional neural networks for predicting wave propagation. In ICLR
Workshop on Deep Neural Models and Differential Equations, 2020.
Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio. Deep learning. MIT press
Cambridge, 2016.
Xiaoxiao Guo, Wei Li, and Francesco Iorio. Convolutional neural networks for steady flow approximation. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge
discovery and data mining, pp. 481–490, 2016.
George Karniadakis and Spencer Sherwin. Spectral/hp element methods for computational fluid
dynamics. Oxford University Press, 2013.
Byungsoo Kim, Vinicius C Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara
Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. In Computer
Graphics Forum. Wiley Online Library, 2019.
Sung-Eun Kim and Ferit Boysan. Application of cfd to environmental flows. Journal of Wind
Engineering and Industrial Aerodynamics, 81(1-3):145–158, 1999.
6

Published as a conference paper at ICLR 2022

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua
Bengio and Yann LeCun (eds.), 3rd International Conference on Learning Representations, ICLR
2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http://arxiv.org/abs/1412.6980.
Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter. Self-normalizing
neural networks. arXiv preprint arXiv:1706.02515, 2017.
Sangseung Lee and Donghyun You. Data-driven prediction of unsteady flow over a circular cylinderusing deep learning. Journal of Fluid Mechanics, 879:217–254, 2019.
Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba. Learning particle
dynamics for manipulating rigid bodies, deformable objects, and fluids. In 7th International
Conference on Learning Representations, ICLR 2019, 2019a.
Yunzhu Li, Jiajun Wu, Jun Yan Zhu, Joshua B. Tenenbaum, Antonio Torralba, and Russ Tedrake.
Propagation networks for model-based control under partial observation. Proceedings - IEEE
International Conference on Robotics and Automation, 2019-May:1205–1211, 2019b. ISSN
10504729. doi: 10.1109/ICRA.2019.8793509.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential
equations. arXiv preprint arXiv:2010.08895, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755–6766, 2020b.
Mario Lino, Chris Cantwell, Stathi Fotiadis, Eduardo Pignatelli, and Anil Anthony Bharath. Simulating surface wave dynamics with convolutional networks. In NeurIPS Workshop on Interpretable Inductive Biases and Physically Structured Learning, 2020.
Wenzhuo Liu, Mouadh Yagoubi, and Marc Schoenauer. Multi-resolution graph neural networks
for pde approximation. In International Conference on Artificial Neural Networks, pp. 151–163.
Springer, 2021.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li Fei-Fei, Joshua B. Tenenbaum, andDaniel L.K. Yamins. Flexible neural representation for physics prediction. Advances in Neural
Information Processing Systems, 2018-Decem:8799–8810, 2018a. ISSN 10495258.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li F Fei-Fei, Josh Tenenbaum, andDaniel L Yamins. Flexible neural representation for physics prediction. In Advances in neural
information processing systems, pp. 8799–8810, 2018b.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W Battaglia. Learning meshbased simulation with graph networks. In (accepted) 38th International Conference on Machine
Learning, ICML 2021, 2021.
Hillel Rubin. Environmental fluid mechanics. CRC Press, 2001.
Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller,
Raia Hadsell, and Peter Battaglia. Graph networks as learnable physics engines for inference andcontrol. 35th International Conference on Machine Learning, ICML 2018, 10:7097–7117, 2018.
Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graphnetworks with ode integrators. arXiv preprint arXiv:1909.12790, 2019.
Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter
Battaglia. Learning to simulate complex physics with graph networks. In International Conference on Machine Learning, pp. 8459–8468. PMLR, 2020.
Nils Thuerey, Konstantin Weissenow, Lukas Prantl, and Xiangyu Hu. Deep Learning Methods for
Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows. AIAA Journal, 58(1):15–26, oct
2018. ISSN 0001-1452. doi: 10.2514/1.j058291. URL http://arxiv.org/abs/1810.
08217.
7

Published as a conference paper at ICLR 2022

Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Accelerating eulerian
fluid simulation with convolutional networks. In International Conference on Machine Learning,pp. 3424–3433. PMLR, 2017.
Steffen Wiewel, Moritz Becher, and Nils Thuerey. Latent space physics: Towards learning the
temporal evolution of fluid flow. In Computer Graphics Forum. Wiley Online Library, 2019.
Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip. A
comprehensive survey on graph neural networks. IEEE Transactions on Neural Networks andLearning Systems, 2020.

ADATASETS DETAILS

A.1

A DVECTION DATASETS

We solved the two-dimensional advection equation using Nektar++, an spectral/hp element solver(Karniadakis & Sherwin, 2013; Cantwell et al., 2015). As initial condition, φ0 , we take an scalar
field derived from a two-dimensional Fourier series with M × N random coefficients, specifically
φ0 =

N
M X
X



cm,n φm,n · exp − 2(x − xc )2 − 2(y − yc )2 ,(3)

m=0 n=0with φm,n = ℜ exp i2π(mx + ny) .
(4)Coefficients cm,n are sampled from a uniform distribution between 0 and 1, and integers M and N
are randomly selected between 3 and 8. In equation (4), xc and yc are the coordinates of the centreof the domain. The initial field φ0 is scaled to have a maximum equal to 1 and a minimum equal to−1. We created training and testing datasets containing advection simulations with 50 time-points
each, equispaced by a time-step size dt = 0.03. A summary of these datasets can be found in Table
3.
Training datasets. We generated two training datasets: AdvBox with 1500 simulations and
AdvInBox with 3000 simulations. In these datasets we impose a uniform velocity fields withrandom values for u and v, but constrained to u2 + v 2 ≤ 1. In dataset AdvBox the domain isa square (x ∈ [0, 1] × [0, 1]) with periodicity in x and y. In dataset AdvInBox the domain is arectangle (x ∈ [0, 1] × [0, 0.5]) with periodicity in y, a Dirichlet condition on the left boundary and ahomogeneous Neumann condition on the right boundary – as an additional constraint, u ≥ 0. During training, a new set of nodes V 1 is selected at the beginning of every iteration. The node count
was varied smoothly across the different regions of the domains. The sets of node were created withGmsh, a finite-element mesher. The element
was set to 0.012 in the corners and the√ size parameter√
centre of the training domains, and set to 10 or 1/ 10 times that value at one random control pointon each boundary. The mean number of nodes in |V1 | for AdvBox and AdvInBox predictions are
9802 and 5009 respectively.
Testing datasets. We generated eight testing datasets, each of them containing 200 simulations.
These datasets consider advection on more complex open and closed domains with non-uniform
velocity fields. The domains employed are represented in Figure 5, and the testing datasets arelisted in Table 3. The velocity fields were obtained from the steady incompressible Navier-Stokes(N-S) equations with Re = 1. In dataset AdvTaylor the inner and outer walls spin at a velocity
randomly sampled between −1 and 1. In datasets AdvCircle, AdvSquare, AdvEllipseH,
AdvEllipseV and AdvSplines there is periodicity along x and y, and a horizontal flow rate
between 0.2 and 0.75 is imposed. The obstacles inside the domains on the AdvSplines dataset are
made of closed spline curves defined from six random points. Dataset AdvCircleAng is similar
to AdvCircle, but the flow rate forms an angle between −45 deg and 45 deg with the x axis. The
domain in dataset AdvInCir has periodicity along y, a Dirichlet condition on the left boundary(with 0.2 ≤ u2 + v 2 ≤ 0.75 and −45 deg ≤ arctan(v/u) ≤ 45 deg), and a homogeneous Neumanncondition on the right boundary. The set of nodes V 1 were generated using Gmsh with an element
size equal to 0.005 on the walls of the obstacles and 0.01 on the remaining boudaries.
8

Published as a conference paper at ICLR 2022

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Physical domains (black areas) on our testing datasets.
Table 3: Advection training and testing datasets
Dataset
AdvBoxAdvInBox
AdvTaylor
AdvCircle
AdvCircleAng
AdvSquare
AdvEllipseH
AdvEllipseV
AdvSplines
AdvIncir

A.2

Flow type

Domain#NodesTrain/Test

Open, periodic in x and yOpen, periodic in y
Closed, Taylor-Couette flow
Open, periodic in x and yOpen, periodic in x and yOpen, periodic in x and yOpen, periodic in x and yOpen, periodic in x and yOpen, periodic in x and yOpen, periodic in y[0, 1] × [0, 1][0, 1] × [0, 0.5]
Figure 5a
Figure 5b
Figure 5b
Figure 5c
Figure 5d
Figure 5eFigure 5f
Figure 5b

9601-10003
4894-5135
7207
19862
19862
19956
20210
20221
19316-20389
19862

Training
Training
Testing
Testing
Testing
Testing
Testing
Testing
Testing
Testing

I NCOMPRESSIBLE FLUID DYNAMICS DATASETS

We solved the two-dimensional incompressible Navier-Stokes equation using the high-order solver
Nektar++. We consider the flow around an infinite vertical array of circular cylinders, with diameter
D = 1, equispaced a distance H randomly sampled between 4D and 6D. The width of the domain
is 7D and the cylinders axis is at 1.5D from the left boundary. The left boundary is an inlet withu = 1, v = 0 and ∂p/∂x = 0; the right boundary is an outlet with ∂u/∂x = 0, ∂v/∂x = 0 and
p = 0; and, the cylinder walls have a no-slip condition. In our simulations, we select Re values thatyield solutions in the laminar vortex-shedding regime, and we only include the periodic stage to ourdatasets. The sets of nodes V 1 employed for each simulation were created using Gmsh placing morenodes around the cylinders walls. The mean number of nodes in these sets is 7143. Each simulation
contains 100 time-points equispaced by a time-step size dt = 0.1. The training and testing datasetsare listed in Table 4.

B

M ODEL DETAILSHyper-parameters choice. The number of edges going to each node was set to k = 6, and thenumber of layers in each MLP (encoder, decoder and edge and node update functions) was three,with 128 neurons per hidden layer. All MLPs (except the decoder) use SELU activation functions(Klambauer et al., 2017), and, batch normalisation (Ba et al., 2016). The grid sizes used for gener-

Table 4: Incompressible flow training and testing datasets Table 5: Cell sizes for coarsening to levels 2, 3
and 4
DatasetRe#Simulations

Train/TestNS
NSMidReNSLowRe
NSMidRe

500-1000
500-1000
100-500
1000-1500

1000
250
250
250

Training
Testing
Testing
Testing

9

Cell size

Advection

Fluid dynamics

d2x , d2y
d3x , d3y
d4x , d4y

0.02
0.04
0.08

0.15
0.30
0.60

Published as a conference paper at ICLR 2022

ating the coarser graphs are collected in Table 5. The number of MP layers we used at each scale
are listed in the Table 6.
Table 6: Number of MP layers at each scale for L = 1, 2, 3 and 4
L

Advection

Fluid dynamics

L=1
L=2
L=3
L=4

M1 = 4
M1 = 2, M2 = 4
M1 = 2, M2 = 2, M3 = 4
M1 = 2, M2 = 2, M3 = 2, M4 = 4

M1 = 8
M1 = 4, M2 = 4
M1 = 4, M2 = 2, M3 = 4
M1 = 4, M2 = 2, M3 = 2, M4 = 4

Training details. We trained MultiScaleGNN models on a internal cluster using 4 CPUs, 86GB ofmemory, and a RTX6000 GPU with 24GB. We fed 8 graphs per batch. First, each training iteration
predicted a single time-point, and, every time the training loss decreased below a threshold (0.01 foradvection and 0.005 for fluid dynamics), we increased the number of iterative time-steps by one, up
to a limit of 10. We used the loss function given by
L = MSE û(t, x1V 1 ), u(t, x1V 1 )+ λd MAE û(t, x1V 1 ∈ ∂D D), u(t, x1V 1 ∈ ∂D D)û(t, x1rk ) − û(t, x1sk ) u(t, x1rk ) − u(t, x1sk )
λe X
+ 1
MSE
,
.
(5)
|E |
||ek ||2||ek ||21∀k∈Ewith λd = 0.25, and, λe = 0.5 for advection and λe = 0 for fluid dynamics. The initial time-point
was randomly selected for each prediction, and, we added to the initial field noise following a uniform distribution between -0.01 and 0.01. After each time-step, the models’ weights were updated
using the Adam optimiser with its standard parameters Kingma & Ba (2015). The learning rate
was set to 10−4 and multiplied by 0.5 when the training loss did not decrease after six consecutive
epochs, also, we applied gradient clipping to keep the Frobenius norm of the weights’ gradients
below or equal to one.

CC OARSENING COMPARISON

An important question in the design of MultiScaleGNN was how to coarsen Vl to Vl+1 . In thefield of numerical simulations several coarsening algorithms have been developed to guarantee thestability and accuracy of the simulations. Liu et al. (2021) and Belbute-Peres et al. (2020) usedcoarsening techniques from numerical simulations and then interpolated the nodes attributes into
the coarsened set of nodes. We found that our cell-grid coarsening and learnt MP from the highto the low-resolution graph (see section 3.2) performs significantly better than such coarseninginterpolation approach. Table 7 shows the MAE on the incompressible flow testing datasets forMultiScaleGNN models using both pooling strategies and a random coarsening, which randomly
drops nodes to match the same |Vl |/|Vl+1 | ratio. In contrast to coarsening algorithms adopted fromnumerical simulations, cell-grid coarsening does not preserve the spatial node density from the original discretisation, but it results into a pseudo-structured discretisation. This discretisation may help
to spread the node information more efficiently on all directions and regions and to ease the learningof the parameters of the MP layers at low resolution scales. Also, instead of a vanilla interpolation,in the DownMP the node information is passed to the low-resolution graph using a learnt MP layer,
and the edge information is also passed to the edges of the low-resolution graph.

10

Published as a conference paper at ICLR 2022

Table 7: MAE ×10−2 on the N-S testing datasets for L = 1, 2, 3, 4 with three different coarsening algorithms
Datasets

L

Cell grid

Guillard’s

RandomStep 99

All

Step 99

All

Step 99

All

NSMidRe

2
3
4

4.759
3.851
3.456

3.663
3.081
2.825

6.756
5.735
4.494

4.393
3.952
3.330

7.191
6.431
6.337

4.741
4.606
4.428

NSLowRe2
3
4

12.34611.707
7.338

8.211
8.203
5.532

9.053
9.099
14.358

6.378
6.946
9.164

10.7368.977
10.335

7.255
6.672
7.521

NSHighRe

2
3
4

7.879
6.980
5.826

9.002
7.096
5.871

10.776
10.331
11.222

7.608
7.556
7.837

11.024
9.631
9.572

7.842
7.28
7.202

11