Full steam ahead!

Now that we understand the physical laws governing fluids — the Navier-Stokes Equations (NSEs) – from Flow #1 and how they can be modelled with Smoothed Particle Hydrodynamics (SPH) from Flow #2, it’s time to focus on how to handle the generated simulation data.

In this post, we explain how to construct the input/output objects necessary to then deploy a Deep Learning model to learn fluid dynamics. This will take us to graphs and their matrix representations.

Once their basics are covered, we construct the graphs, with appropriately selected features from the SPH data. This will get us closer to generating a dataset from which the model can be trained.

So, let’s flow!

Flow #3:
Representing Smoothed Particle Hydrodynamics (SPH) Data with Graphs

The world
… of graphs

In Particle Land, SPH particles are nodal points where quantities of interest are evaluated. In Flow #2, we have seen how any function \(f\) is evaluated locally is made by averaging information over neighbours.

Unsurprisingly, a very suitable object that can naturally encapsulate this information is a graph. Drily put, a graph is a data structure composed of a set of nodes and a set of edges that connect pairs of nodes [5].

Graphs can model almost any relationship. Nodes usually represent entities of interest in the problem domain — such as people, cities or molecules. Edges encode connections or relationships between the entities, such as two people being acquainted, cities connected by a road or chemical reactions between molecules [1].

The edges can be made directional too. If node \(i\) and \(j\) are connected both ways, the edge does not have a specified direction. In this situation, the graph is undirected. Otherwise, if node \(i\) is connected to node \(j\), but not vice-versa, then the graph is directed, as the graph in Fig. 1:

Fig. 1: A directed graph. Node 1 is connected to nodes 3 and 4, node 3 is connected to node 4 and nodes 2 and 4 are connected to each other. Credits: David Carvalho / Inductiva.

But how can we encode the structure of a graph mathematically?

Representing
… graphs

The connectivity structure of a graph containing \(n\) nodes may be represented by an \(n \times n\) matrix, the adjacency matrix \(A\). Its entries \(A_{ij}\) indicate whether there is a connection from node \(i\) to node \(j\) (\(A_{ij} = 1\)) or not (\(A_{ij} = 0\)) [2].

For example, the adjacency matrix \(A\) of the graph in Fig. 1 is:

\[A = \begin{pmatrix} 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \end{pmatrix}\]

If the graph is undirected, the adjacency matrix is naturally symmetric, since all connections go both ways.

Featuring
…features

Fine, nodes and edges — but so far, no information is present. This can be included in three levels:

• the state: the graph can have \(u\) global features pertaining to the entire system — think of the average score amongst all pupils or the centre of mass of all particles. They are represented by a global feature vector \(\mathbf{u} \in \mathbb{R}^u\);
• the entities: At the node level, several features can be associated — think of a person’s age or height, a city’s population (… or the particle velocity in a fluid). The \(x\) features of node \(i\) are represented by a feature vector \(\mathbf{x}_i \in \mathbb{R}^{x}\);
• the relations: Much like the entities, relationships between them can also have features — think of the distance between cities (…or particles!). The \(e\) features of an edge connecting node \(i\) to node \(j\) are represented by a feature vector \(\mathbf{e}_{ij} \in \mathbb{R}^{e}\).

Fig. 2: Graph with the associated feature vector structures at the global (\(\mathbf{u}\)), edge (\(\mathbf{e}_{ij}\)) and node (\(\mathbf{x}_j\)) level. Credits: David Carvalho / Inductiva.

These vector objects can be stored in matrices. The features of all nodes in a graph can then be aggregated in a feature matrix \(X\) of size \(n \times x\), whose \(i^{\rm th}\) row is the feature vector \(\mathbf{x}_i\) of node \(i\).

Likewise, the edge features of all \(n_e\) edges in a graph can be aggregated in the edge feature matrix \(E\), of size \(n_e \times e\).

To fix an order of \(E\), its rows can be sorted in the Coordinate Format (COO) [3]. In it, edges are first sorted by row index and then by column index of the entries, as set by the in the adjacency matrix \(A\).

Constructing graphs
… for SPH data

Fine — graphs can store a lot of information about a complex system — but what motivates the use of graph data structures in this problem?

Since nodes can be linked in a custom way and carrying custom information, learning from graphs introduces a strong relational bias between the items, which would not be possible if simple, tabular data was used instead.

As we saw in the previous post, the dynamics of a particle is not solely determined by its own properties, but also by the ones from neighbouring particles.

Well, it makes sense to use a graph whose nodes represent fluid particles and edges represent physical interactions between them, including the notion of neighborhood.

Fig. 3: The position of the particle nodes in various instants. Here, a column of water is thrown against a wall — Splash!. Credits: Pedro Pereira / Inductiva.

Outstanding — particle data can be encoded in graphs.
But what do we want to predict exactly?

Finding
… a target output

In order to apply supervised Machine Learning to graph structures, we must specify what our target output \(Y\) is for a given input.

The graph input is then given by these matrix representations — \(A\), \(X\) or \(E\). For graphs, three distinct kinds of tasks exist depending on the type of task.
For a number of features \(y\), and general output \(Y\) of the task can happen at the:

• graph-level prediction, where a global feature vector \(Y \in \mathbb{R}^y\) can be learnt
• node-level prediction, where labels for each node can be stored in a matrix \(Y \in \mathbb{R}^{n \times y}\)
• edge-level prediction, where each edge can be stored in a matrix of size \(Y \in \mathbb{R}^{n_e \times y}\).

In this situation, we are interested in building a Deep Learning (DL) model that, given an input graph, can understand the underlying physics governing particle trajectories and determine the time (or temporal) evolution of the system.

The particle acceleration \(\mathbf{a}_j\) is a good target. If \(\mathbf{a}_j\) is known in a given timeframe, their velocities \(\mathbf{v}_j\) and positions \(\mathbf{r}_j\) can be computed in time by an integration scheme. The consequence? Given enough steps, the full trajectories throughout each timeframecan be computed!

So the Input/Output we prescribe is as follows: given a graph \(\mathcal{G}^t = (A^t, X^t, E^t)\) representing the system at time \(t\), we are interested in learning the particle acceleration vector \(\mathbf{a}^t\), i.e. the matrix \(Y = \mathbf{a}^t \in \mathbb{R}^{n \times 3}\) whose i ^th row stores the three components of the acceleration felt by particle in node \(i\).
Note that we will not need to assign any global features \(\mathbf{u}\).

The mapping seems clear now. But which features do we use for the graph inputs?

Finding
… the right features

Choosing the right features is essential to ensure a model can learn (even for the ones used in Deep Learning!) — so let us go step-by-step.

Throughout the Flow series, we follow closely the model proposed by Sanchez-Gonzalez et al [8] — a fantastic read!

After a SPH simulation is run [4], the data generated comprises the particle positions \(\mathbf{x}^t\) for each time \(t\) up to a total simulation time \(T\). Each timeframe contains data which can be converted to a graph. Let’s pick a certain timeframe \(\mathbf{x}^t\).

… For the edges

Before we define any features, we must compute which edges are present. We are interested in the neighborhood of each node, meaning that edges are added between nodes \(i\) and \(j\) if their distance is less than a specified connectivity radius \(R\), i.e.

\[||\mathbf{r}^{t}_i − \mathbf{r}^{t}_j || \leq R\]

The connectivity radius \(R\) is analogous to the SPH kernel smoothing length \(h\) that we discussed in Flow #2. This graph operation may be performed by standard nearest neighbour search algorithms, such as the kd-tree algorithm [7].

Phew — the adjacency matrix \(A^t\) has been constructed.
We now need to assign edge feature vectors \(\mathbf{e}^t_{ij}\). We use the relative displacement vector and corresponding distance:

\[\mathbf{e}^t_{ij} = (\mathbf{r}^{t}_i − \mathbf{r}^{t}_j, ||\mathbf{r}^{t}_i − \mathbf{r}^{t}_j||)\]

Woohoo — the edge feature matrix \(E^t\) is also constructed.

… For the nodes

At the node level, we need to include enough previous dynamical information so the prediction is successful. We do this by featuring the previous \(C\) velocities i.e. the sequence \(v^{t-C + 1}_j, ..., v^{t}_j\), which can be calculated using a Finite Difference estimate:

\[\begin{align} \mathbf{v}^t \approx (\mathbf{x}^t − \mathbf{x}^{t-1})/ \Delta t && \mathbf{a}^t \approx (\mathbf{x}^{t+1} − \mathbf{x}^{t})/ \Delta t^2 \end{align}\]

\(C\) will be a hyperparameter controlling how much past information is needed.

Particles close to the boundary undergo more complex dynamics e.g. due to backscattering or corners. The information about the simulation box is encoded in a 6-dimensional vector \(\mathbf{d}_i\) whose components store the the particle distance to each of the six cubic walls of the domain. The node representing particle \(i\) thus has as the feature vector:

\[\mathbf{x}^{t}_i = (\mathbf{v}^{t−C+1}_i, \dots , \mathbf{v}^{t}_i , \mathbf{d}^t_i)\]

And we finally have the (node) feature matrix \(X^t\) assembled. Our graph \(\mathcal{G}^t\) is ready to serve!

A graph is constructed for each time frame of a simulation — so there can be quite a lot of them to handle! The number of timeframes in a simulation is determined by \(T\) and the time step \(\Delta t\).

For instance, a fluid flow simulated for \(T = 4\) s with \(\Delta t = 5\) ms results in \(T/\Delta t = 800\) timeframes! If, say, \(C = 5\) previous velocities are used, a total of \(800 = 796\) graphs can be built with a past window spanning at least 5 timeframes.

In the
next post…

We now understand how to encapsulate the flow data obtained from a SPH simulation in a graph \(\mathcal{G}^t\). The task has also been fixed: we want to predict \(\mathbf{a}^t\).

But how exactly can we compute the desired mapping \(f : \mathcal{G}^t \mapsto \mathbf{a}^t\) ?
Needless to say, this is exactly what our Deep Learning model will be responsible for.

In the next post, we will dive deeper into the model architecture and design choices that will ultimately will allow us to compare ground truth (set by SPH data) and our own very model!

Until then…
~~~ flow away! ~~~

---

References
… & Remarks

[1] A hands-on introduction to graphs with a lot of concise examples.

[2] A collection of many types of graphs (from simple graphs to digraphs and networks) with their corresponding adjacency matrix \(A\).

[3] COO is particularly useful in constructing sparse matrices in a memory-efficient way. Here are some methods in Python’s SciPy.

[4] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia, “Learning to simulate complex physics with graph networks” in International Conference on Machine Learning. PMLR, 2020, pp. 8459–8468.

[5] Reducible — A pretty neat introduction to Graph Theory (from a Computer Science perspective).

[6] The Flow series closely follows the great work that was done for the Masters thesis by Pedro Pereira, one of our student collaborators. You can access it in the link!

[7] S. Skiena, The Algorithm Design Manual, ser. Texts in Computer Science. Springer International Publishing, 2020, ch. 15.6, 20.5.

[8] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia, “Learning to simulate complex physics with graph networks” in International Conference on Machine Learning. PMLR, 2020, pp. 8459–8468.