#
Flow #3:

Representing Smoothed Particle Hydrodynamics (SPH) Data with Graphs

**Pedro Ferro Pereira**Co-author

**David Carvalho**Co-author | Editor

**Fábio Cruz**Technical Reviewer

**Ivan Pombo**Reviewer

*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 **I**nput/**O**utput 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.

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:

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:

\(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.

## Recent posts from our blog

Sofia Guerreiro

Cristiana Carpinteiro

In this series of blog posts we will explore a specific case of the use of AI in the pharmaceutical industry - using Graph Neural Networks for predicting binding affinity. But for now, let’s start by understanding the problem of drug discovery and some fundamental concepts like binding affinity.

Ivan Pombo

Luís Sarmento

Coastal engineering projects protect the coast from erosion, flooding and other events. Due to their cost, the design phase of these projects is heavily based on computational simulations. In this blog post, we enlighten how Inductiva API allowed researchers to scale their coastal engineering simulations.