# A fluid simulation system based on the MPS method

André Luiz Buarque Vieira e Silva
Voxar Labs - Centro de Informática
Universidade Federal de Pernambuco
Recife/PE 50740-560, Brazil.
<EMAIL_ADDRESS>
&Caio José dos Santos Brito
Voxar Labs - Centro de Informática
Universidade Federal de Pernambuco
Recife/PE 50740-560, Brazil.
<EMAIL_ADDRESS>
&Francisco Paulo Magalhães Simões
Voxar Labs - Centro de Informática
Universidade Federal de Pernambuco
Recife/PE 50740-560, Brazil.
Departamento de Informática
Instituto Federal de Pernambuco, Campus Belo Jardim
Belo Jardim/PE 55145-065, Brazil
<EMAIL_ADDRESS>
&Veronica Teichrieb
Voxar Labs - Centro de Informática
Universidade Federal de Pernambuco
Recife/PE 50740-560, Brazil.
<EMAIL_ADDRESS>

###### Abstract

Fluid flow simulation is a highly active area with applications in a wide
range of engineering problems and interactive systems. Meshless methods like
the Moving Particle Semi-implicit (MPS) are a great alternative to deal
efficiently with large deformations and free-surface flow. However, mesh-based
approaches can achieve higher numerical precision than particle-based
techniques with a performance cost. This paper presents a numerically stable
and parallelized system that benefits from advances in the literature and
parallel computing to obtain an adaptable MPS method. The proposed technique
can simulate liquids using different approaches, such as two ways to calculate
the particles’ pressure, turbulent flow, and multiphase interaction. The
method is evaluated under traditional tests cases presenting comparable
results to recent techniques. This work integrates the previously mentioned
advances into a single solution, which can switch on improvements, such as
better momentum conservation and less spurious pressure oscillations, through
a graphical interface. The code is entirely open-source under the GPLv3 free
software license. The GPU-accelerated code reached speedups ranging from 3 to
43 times, depending on the total number of particles. The simulation runs at
one fps for a case with approximately $200,000$ particles. Code:
https://github.com/andreluizbvs/VoxarMPS

_K_ eywords MPS $\cdot$ Framework $\cdot$ Numerical improvements $\cdot$ Fluid
models $\cdot$ Parallelization

PROGRAM SUMMARY

Program Title: Voxar MPS
Licensing provisions: GNU General Public License, version 3
Programming language: C++ and CUDA
Computer: Tested on CPUs: Intel Core i7-6820HK and Intel Core i7-7820HK and
GPUs: GTX 1080 (mobile) and GTX 1080 Ti
Operating System: Windows 10
CUDA: Tested on version 10.1 with driver version 445.87
Has the code been vectorised or parallelised: Different threads of CPU or
number of cores of GPU
RAM: Tens of MB to several GB, depending on the scenario
Nature of problem: The Voxar MPS code has been developed to study the flow of
incompressible fluids that requires high computational cost.
Solution method: Voxar MPS is an implementation of the Moving Particle Semi-
implicit, a Lagrangian meshless particle method for incompressible fluids.

##  1 Introduction

Some of the most common problems in naval hydrodynamics involve the study of
fluid flow. For this, it is necessary to deal with large deformations near the
surface, such as those presented in a good portion of computational mechanics
problems [1]. Also, Computer Graphics (CG) and Virtual Reality (VR)
applications are constantly turning more realistic and interactive/real-time
fluid simulations have been a frequent research topic in that area [2, 3, 4].
Conventional techniques, as the Finite Element Methods (FEM) and Finite
Difference Methods (FDM) are relatively inefficient when dealing with large
deformations [5] [6]. As an alternative, there are the Lagrangian meshfree or
the particle-based methods. They achieve flexibility in situations where the
classic techniques are too complex [7].

Two of the meshless techniques are the Moving Particle Semi-implicit method
(MPS) [8] and the Smoothed Particle Hydrodynamics (SPH) [9] [10], the latter
initially intended to astrophysics applications and then adapted to fluid
simulation. The MPS authors idealized it to simulate the flows of
incompressible fluids, which refers to a fluid whose material density is
constant within a fluid parcel, a property found in liquids. Its main
difference from the original SPH, which can be considered an advantage for the
MPS regarding numerical precision, is that the calculations adopt a semi-
implicit predictor-corrector model [8]. However, the SPH is more prevalent in
CG and VR applications due to the high computational load occasioned by the
MPS calculations, including solving the Poisson Pressure Equation (PPE).

In [11], the authors highlight the current achievements and future
perspectives for projection-based particle methods, which are the ones that
require solving a PPE. In it, the authors present a set of papers regarding
the applicability of the MPS in ocean engineering, including wave breaking
[12] [13], wave overtopping [14], wave impact [15] [16], green water on ships
[17], sediment transport [18], waves generated through landslide [19] and
fluid-structure interactions [20] [21].

One of the main issues of meshfree methods, in general, are in the spurious
pressure oscillation of the particles [15]. Another disadvantage of particle-
based approaches that rely on numerical precision is the application runtime
[22]. Works found in the literature provide various kinds of performance
speedups using a General Purpose Graphics Processing Unit (GPGPU) but always
focusing on the standard MPS [23, 22, 24].

The system developed in this work allows combining numerical improvements to
the MPS with different types of fluid flows. These include multiphase and
viscoplastic simulations and parallel computation at different levels:
multithreaded CPU and GPU. A fine-tuning of these combinations are available
through a graphical user interface. It also allows the simulation of dam break
cases, flood simulations, oil spilling disasters, and others with decent
precision in less time. These features help to achieve general requirements in
which precision and or time are somewhat relevant factors.

This research aims to provide a numerically enhanced MPS method implementation
compared to its standard one, as well as expand the simulation tuning
possibilities for different types of fluids and situations. Additionally, CPU
and GPU parallelization are also provided, so more significant cases can reach
interactive simulation and real-time rates. Modify and add new models to the
system is straightforward. The resulting solution can assist in engineering
problems regarding natural and environmental disasters in coastal and
floodable areas. It can also be used in CG and VR, and, for some cases, real-
time applications.

An unprecedented parallelization of an improved MPS version coupled with fluid
models, such as turbulent and multiphase flow, is developed to run in CPU,
through OpenMP [25], and in GPU, by using CUDA [26]. The tunable parameters
are:

  1. 1.

Multiple computational methods such as different pressure calculation
approaches, different values for the kinematic viscosity, turbulent flow, and
multiphase interaction;

  2. 2.

Improvements such as better momentum conservation, more precise discretization
of differential operators and less spurious pressure oscillations;

  3. 3.

Choose between parallelized implementations using OpenMP and CUDA;

  4. 4.

Interactive rates for simulations containing $10^{5}$ particles.

The code is open-source, available to the community under the GNU General
Public License v3.0 (GPLv3) license. There is also a README file with a
detailed step-by-step installation description and information on how to use
the developed system.

In the second section, there is a contextualization on the area, presenting
modifications to the method over the years and acceleration techniques. Then,
there is a section explaining the used technique and presenting a set of
variations, improvements, and applications. Afterwards, there is an
explanation regarding implementation, acceleration structures and user
interface. Next, the main results are showcased, presenting each test case
utilized, and their purpose. In the test cases section, a discussion is
performed through a comparison with the results found in the literature.
Subsequently, there is a performance analysis of CPU and GPU versions, memory
usage, the speedups provided by the GPU parallelization, and the simulation
frame rate. Finally, there is a discussion of final remarks and future works.

##  2 Related works

Through the years, disasters involving natural phenomena have triggered
research in different areas on how to avoid them. Fluid simulation focusing on
liquids is one of these areas. To be able to simulate liquids coherently,
meshless methods usually model the fluid flow as Weakly Compressible (WC) and
Fully Incompressible (FI), which guarantees less fluid density oscillations.
The MPS has also been providing significant assistance in that field since it
intrinsically simulates incompressible flows [27]. The textbook chapter of
Almeida et al. [28] is an introduction on the subject since it provides a
thorough explanation of the MPS as well as some variations and implementation
details.

Despite the high level of applicability of the MPS, as shown in the previous
section, there are only a few solutions or software that use this method, such
as MPS-RYUJIN [29], Particleworks [30] and MPARS [31]. Refs. [29] [30] are not
open-source and Ref. [31] does not provide any kind of parallelization, as
opposed to the present work. Although there are other frameworks for fluid
simulation, some of them being widely used, such as the DualSPHysics [32] and
GPUSPH [33], they mainly focus on the SPH and its variations.

###  2.1 Modifications to the method

Similar to the other meshfree methods, the MPS technique suffers from
instability problems. Some of these issues are related to numerical errors at
the boundaries, i.e., free-surfaces or solid boundaries interactions. Some
works document and discuss the reasons for those problems [34] [16]. As an
attempt to overcome these issues, some authors proposed changes to the method
throughout the years; one of the most critical issues with the MPS is the
spurious pressure oscillation. This subsection presents some of the works that
successfully addressed those issues.

A set of papers by Khayyer & Gotoh presents valuable insights and improvements
to this problem, most of them proposing corrected differential operator models
(Laplacian and gradient). [15] proposes modifications to the MPS to diminish
spurious pressure fluctuations. The authors introduce a new formulation of the
source term of the PPE, which is referred to as a Higher order Source term
(HS), thus creating the CMPS-HS after combining this modification with their
previous work. The authors in [35] focus on the Laplacian model used in the
MPS. They derive a Higher order Laplacian model (HL) for the discretization of
the Laplacian operator to refine further and stabilize the pressure
calculation. Both Laplacian of pressure and the one corresponding to the
viscous forces benefit from this model. The authors remarked that, although
the modifications improve pressure calculations, the numerical results still
presented some unphysical numerical oscillation during tests.

Khayyer & Gotoh [36] present two new modifications to resolve the shortcomings
present in the method from their previous work. The first improvement deals
with unphysical numerical oscillation caused by the source term in the PPE.
Due to that issue, additional terms were added to it, referred by Error
Compensating parts in the Source term of the PPE (ECS). The other change deals
with situations with tensile instability [37]. It consists of a corrective
matrix inserted in the pressure gradient calculations to achieve a more
accurate approximation of the differential operator in question. The
modifications here presented are going to be detailed further in section 3.

In [38], the authors show the latest advances related to particle methods
applied to coastal and ocean engineering, where they also include fully
explicit methods. In one of them, Tayebi & Jin [39] proposes the Moving
Particle Explicit (MPE) that solves an equation of state in a fully explicit
form to obtain the particles’ pressure. Other recent applications of the MPS
brought up by the referred review are related to oil spilling [40] and the
swash beach process [41]. In the present work, a fully explicit version of the
MPS is studied, implemented, and accelerated.

###  2.2 Acceleration

Since the MPS is fully meshless, the particles are not connected explicitly by
any edge, so, it is possible to optimize some computational aspects of the
simulation, such as by parallelization, by cluster technology or General
Purpose GPU (GPGPU) techniques [42]. Shakibaeinia et al. [43] [44] presented
an alternative form of the MPS, where they incorporate more straightforward
and effective formulations and fluid behaviours/flows models into the method.
One of the advantages of doing so is to benefit from the MPS’ numerically
precise formulations without the main efficiency bottleneck, which is the
Poisson pressure equation’s solution.

Shakibaeinia and Jin [43] developed a technique based on the MPS interaction
model in which a WC model replaces the FI one. By doing this, they decrease
the necessary computation time. In Ref. [44], the authors extended their
previous work by proposing a straightforward model of an immiscible multiphase
method. They also employ and investigate different techniques for the
viscosity model. In order to address turbulence issues in wave dynamics, they
use the Large Eddy Simulation (LES) concept to formulate a sub-particle scale
(SPS) turbulence model. The modifications mentioned above are studied and
taken into consideration in the present work.

Hori et al. [23] developed a GPU-accelerated version of MPS using CUDA. The
authors focused on the search of neighbouring particles and the iterative
solution of the linear system generated by the PPE, which generates a
considerable computational load. They optimize the neighbouring particles
search through a cell grid, in which there is a specific cell for each
particle according to the particle’s position. To compare accuracy and
performance between the CPU and GPU-based codes, they execute 2-dimensional
calculations of an elliptical drop evolution and a dam break flow. Finally,
the reported speedup achieved is about 3 to 7 times. These speedup rates serve
as comparison reference here.

In Ref. [24], the authors focus on reducing the MPS runtime by replacing how
the method calculates the particles’ pressure. Instead of using a PPE, they
implement an equation of state, similar to the work of Shakibaeinia and Jin
[43]. It also aims to accelerate it by exploring its parallelization potential
through a multi-core CPU, single-node GPU, and a multi-node GPU cluster
environment. The authors use a domain subdivision approach to enable a
simulation with a higher number of particles. For a 3D dam break test with
$~{}700,000$ particles, the OpenMP solution could reach 5.3 times speedup,
while the single-node GPU could reach 14.5 times speedup compared to a single-
threaded CPU execution. The multi-node GPU with nine processes performs
approximately 5.5 times faster than the single-node GPU. The authors claim
that the proposed algorithm allows extensive WC-MPS simulations in distributed
memory systems with reduced communication overhead. These speedup values are
also a comparison reference to this work.

Unlike the works here presented, this work not only focuses on applying
parallelization to the MPS calculations but also applying it to the added
improvements and modifications. This combination of features is unprecedented.
Besides, both parallelized versions (OpenMP and CUDA) can achieve speedups
regarded to the MPS variations developed. Exemplifying, one set up of features
can praise for stability and higher accuracy, while another may seek for
higher performance.

##  3 The Moving Particle Semi-implicit Method

The MPS uses discrete elements called particles in which each of them carries
a set of physical quantities. Since the fluid flow governing equations are for
continuous domains, the continuous differential operators, such as the
derivative, gradient, and Laplacian, need to be discretized. The MPS proposes
discretized models to these operators. In this section, the method is detailed
by showing its governing equations, discretized differential operations, and a
set of variations and improvements to the standard MPS.

###  3.1 Standard method & Governing equations

Koshizuka and Oka [8] models the fluid as a set of interacting particles, in
which their motion is determined through the interaction with neighbouring
particles and the governing equations of fluid motion. Equation 1 and Equation
2 are the continuity equation and Navier-Stokes equation, respectively, which
describe the motion of one viscous fluid flow .

$\frac{1}{\rho}\frac{D\rho{}}{Dt}+\nabla{}\cdot\textbf{u}=0$ (1)

$\frac{\partial\textbf{u}}{\partial
t}=-\frac{1}{\rho}\nabla{p}+\nu\nabla^{2}{\textbf{u}}+\frac{1}{\rho}\nabla\cdot\vec{\tau}+F_{ext}$
(2)

where u is the fluid velocity vector, $t$ is the time, $\rho$ is the fluid
density, $p$ is the pressure, $\nu$ is the laminar kinematic viscosity,
$\boldsymbol{\tau}$ is the sub-particle scale (SPS) or turbulence contributed
by unresolved small motions (detailed in subsubsection 3.2.2) and $F_{ext}$
represent external forces like gravity.

In this method, the domain is discretized into particles, as mentioned above.
They interact with its neighbours through a kernel function $w(r)$, $r$ being
the distance between two particles. A larger kernel size implies in an
interaction with more particles. The system developed in this work provides a
set of kernel functions to choose between them. However, it is recommended
that the kernel function proposed by [45] should be used when running a PPE
(FI version), as in Equation 3.

$w\left(|{}\textbf{r}_{j}-\textbf{r}_{i}|{}\right)=\left\\{\begin{array}[]{l}\frac{r_{e}}{\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}}-1\
\ ,\ \ 0\leq{}\ r<r_{e}\\\ 0,\ \ \ r_{e}\leq{}r\end{array}\right.$ (3)

where $r_{e}$ is the radius of the interaction area and $\textbf{r}_{i}$ and
$\textbf{r}_{j}$ are the positions of particles $i$ and $j$, respectively.
When calculating the pressure through an equation of state (WC version), it is
recommended to use and the kernel function proposed by [43], as in Equation 4.
Those pressure calculation models are addressed further in this section.

$w\left(|{}\textbf{r}_{j}-\textbf{r}_{i}|{}\right)=\left\\{\begin{array}[]{l}\Big{(}1-\frac{\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}}{r_{e}}\Big{)}^{3}\
\ ,\ \ 0\leq{}\ r<r_{e}\\\ 0,\ \ \ r_{e}\leq{}r\end{array}\right.$ (4)

Equation 5 defines $n_{i}$, the particle number density, at the particle’s
position $\textbf{r}_{i}$, which is proportional to the neighbors number of
$i$.

$n_{i}=\sum_{j\not=i}w\left(\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}\right)$
(5)

The continuity equation is satisfied if the particle number density remains
constant, and this constant value is denoted by $n_{0}$.

To identify a free-surface particle, the particle number density of the $i$th
particle just needs to satisfy the condition presented in Equation 6 since on
the free-surface the particle number density drops abruptly.

$n_{i}<\beta{}n_{0}$ (6)

where $\beta$ is constant between $0.8$ and $0.99$. The bigger $\beta$ is, the
bigger will be the number of particles recognized as free-surface. Koshizuka
and Oka [8] recommend to set it to $0.97$, and that is the value adopted here.

Generally, in meshless methods, there are two main approaches to calculate the
particles’ pressure when simulating liquids, the weakly compressible (WC)
approach and the fully incompressible (FI) one, where each one of them has its
advantages and disadvantages. This work implements both approaches, providing
a wide variety of features so the users can decide based on their needs. After
obtaining the particles’ pressure, it is possible to compute the pressure
gradient, enabling the calculations of the velocity values to update the
particles’ positions through a first-order Euler integration.

####  3.1.1 Equation of state

The WC model prioritizes performance since it severely diminishes
computational load in exchange for numerical precision. In the work of [43],
the traditional FI model is replaced by a WC one because assembling and
solving the PPE in each step takes a considerable amount of computation time:
about two-thirds of each time step for a simulation with an order of magnitude
of $10^{3}$ particles. The mentioned work replaces the PPE by an explicit
relation, specifically an equation of state described by [46] and modified by
[47] that is shown below.

$p_{i}^{k+1}=\frac{\rho{}c_{0}^{2}}{\gamma{}}\left({\left(\frac{n_{i}^{*}}{n_{0}}\right)}^{\gamma{}}-1\right)$
(7)

where $p_{i}^{k+1}$ is pressure value of particle $i$ in timestep $k+1$ and
the typical value used for $\gamma=7$. $c_{0}$ is the speed of sound. This
study, in fact, shows a decrease in process time per time step while the
simulation characteristics remains similar to the fully incompressible
approach. Authors refer to this modified MPS as WC-MPS.

This work adopts the approach proposed by [43] due to its similarity to the
weakly compressible approach adopted in the popular WCSPH [48].

####  3.1.2 Poisson pressure equation

As opposed to the weakly compressible approach, the fully incompressible model
calculates the particles’ pressure with higher accuracy even though this leads
to a high computational load. In this case, it is necessary to solve the
Poisson Pressure Equation, shown in Equation 8, which yields a linear system
of equations of the type shown in Equation 9.

$\nabla{}{{}^{2}}p_{i}^{k+1}=-\frac{\rho{}}{\Delta{}t^{2}}\frac{n_{i}^{*}\
-n_{0}}{n_{0}}$ (8) $Ax=b$ (9)

where $A$ is a sparse square matrix of size $N\times N$ which N is the total
particle number in the simulation, the Right-Hand Size (RHS) vector $b$ of
size $N$ stores the source terms and $x$, also of size $N$, represents the
desired pressures of the particles.

The assemble of the coefficient matrix $A$ makes use of the discretized
Laplacian model, shown in Equation 10.

${\nabla{}}^{2}{\varphi{}}_{i}=\frac{2D_{S}}{n_{0}\lambda{}}\sum_{j\not=i}\left({\varphi{}}_{j}-{\varphi{}}_{i}\right)w\left(\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}\right)$
(10)

where $\varphi$ is some physical quantity, $\lambda$ is the weighted average
of the squared distance between particles $i$ and $j$ (or $r_{ij}^{2}$), as
shown in [8].

As previously mentioned in subsection 3.1, this work uses the ICCG to solve
the PPE, like in [8]. It consists of submitting the coefficient matrix $A$
into the incomplete Cholesky factorization, which is generally used as
preconditioning for iterative numerical methods. Afterwards, the conjugate
gradient is applied to solve the linear system of equations iteratively.

###  3.2 Fluid flow models

The improvement of numerical stability and the addition of different models of
fluid flow enable a more reliable simulation in certain situations. This
section presents the used set of models of fluid flows in this work.

####  3.2.1 Multiphase flow

A model that significantly increases the number of possible applications for
the MPS is the one that supports multi-density fluids interaction, the
multiphase flow.

Shakibaeinia and Jin [44] proposed a straightforward multiphase model based on
the MPS, in which it treats the system as a multi-density multi-viscosity
fluid. The model is only applied to a WC-MPS [43], solving a single set of
equations for all phases. In this model, the density differences of particles
of different phases are automatically taken care of, since that, when
calculating a particle’s velocity, its density appears directly in the
equations.

The main issue of this approach arises when dealing with the density
discontinuity near the interface between the fluids, which can result in
pressure field discontinuity. The strategy followed was to use a smoothed
value of density $\langle\rho\rangle_{i}$ instead of the real particles’
density, set for each particle before starting the simulation. The density of
an individual particle is necessary for calculating its pressure. Equation 11
shows the multi-density pressure term in the weakly compressible model applied
in this study.

$\frac{1}{\rho_{i}}\left\langle\nabla
p\right\rangle_{i}=\frac{d}{n_{0}}\sum_{j\neq
i}\bigg{(}\frac{(\rho_{j}/\rho_{i})\alpha_{j}-\alpha_{i}}{r_{ij}^{2}}\mathbf{r}_{ij}w(\left|r_{j}-r_{i}\right|)\bigg{)}$
(11) $\alpha_{i}=\frac{c^{2}_{0}}{\gamma}((\left\langle
n^{*}\right\rangle_{i}/n_{0})^{\gamma}-1)$ (12)

where $d$ refers to the number of dimensions in the simulation and the typical
value of $\gamma=7$ is used, as in Tait’s equation of state [49].

The multiphase model proposed by Shakibaeinia and Jin [44] can only be used by
calculating the particles’ pressure through the equation of state, which is
the WC approach. Despite that, this work adopts this model for that
incompressibility model, given its relative simplicity and stability.

####  3.2.2 Turbulent flow

To calculate the influence of the turbulence term $\boldsymbol{\tau}$,
referred to the unresolved small motion term in [50], the large eddy
simulation (LES) mathematical model for turbulence [51] [52] was employed.
According to the original LES conception, eddies capable of being resolved by
the computational grid are allowed to evolve according to the Navier-Stokes
equations, and a model is employed to represent the turbulence at sub-grid
scales (mesh-based techniques). A sub-particle scale (SPS) model was made
necessary for meshless methods. By introducing the turbulence eddy viscosity
$\nu_{t}$, the unresolved SPS turbulence stress $\tau_{ij}$ in Equation 2 can
be written as shown in Equation 13.

$\frac{\tau_{ij}}{\rho}=2\nu_{t}S_{ij}-\frac{2}{3}k\delta_{ij}$ (13)

where $\delta_{ij}$ is Kronecker’s operator; and $S_{ij}$ is the strain rate
and $k$ is the turbulence kinetic energy, which can be incorporated into the
pressure term when solving the momentum equation Equation 2. The widely used
model by [51] is employed here to formulate the turbulence eddy viscosity.

###  3.3 Numerical improvements

In this section, there is a description of the implemented MPS variations. It
is noteworthy that the universe of MPS variations is vast, and the ones that
were selected stand between improvement impact level and implementation cost.
These variations allowed a version considered sufficiently stable and
physically accurate to be achieved. It is also shown other modifications to
the standard method, which expand the range of applications of the MPS.

####  3.3.1 Momentum conservation

A simple way to achieve consistent conservation of linear momentum is to
ensure a better discretization of the gradient model. Equation 14 shows the
suggested alteration in the pressure gradient formulation by Khayyer and Gotoh
[53].

$\displaystyle\nabla{}\varphi_{i}=\frac{D_{S}}{n_{0}}\Bigg{(}\sum_{j\not=i}\frac{\left({\varphi_{i}+\varphi}_{j})-({\hat{\varphi}}_{i}+{\hat{\varphi}}_{j}\right)}{{\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}}^{2}}\left(\textbf{r}_{j}-\textbf{r}_{i}\right)w\left(\left|{}\textbf{r}_{j}-\textbf{r}_{i}\right|{}\right)\Bigg{)}$
(14)

When the anti-symmetric Equation 14 is applied, linear momentum is conserved.
This method is referred to by the authors as Corrected MPS (CMPS).

####  3.3.2 Pressure calculation

One of the major issues of the MPS, and consequently widely explored, is the
spurious pressure oscillation. Works that presented substantial improvements
in this area, making few and simple modifications to the method, have been
proposed [15, 35]. The authors call the first one the MPS with a Higher order
Source term (MPS-HS) since it presents a new formulation for the calculation
of the derivative of the particle number density $(\frac{Dn}{Dt})$. Using this
variation, the Equation 8 is replaced by the Equation 15 [54].

$\nabla{}{{}^{2}}p_{i}^{k+1}=-\frac{\rho{}}{n_{0}\Delta{}t}\left(\sum_{i\not=j}\frac{r_{e}}{r_{ij}^{3}}\left(x_{ij}u_{ij}\right.\right.+\left.\left.y_{ij}v_{ij}+z_{ij}w_{ij})\Bigg{)}^{*}\right.\right.$
(15)

where $r_{ij}$ is the distance between particles $i$ and $j$. $x_{ij}$,
$y_{ij}$ and $z_{ij}$ represent the distance between particles $i$ and $j$ in
each dimension and $u_{ij}$, $v_{ij}$ and $w_{ij}$ the velocity difference of
particles $i$ and $j$ in each dimension.

Another improvement to the implemented pressure calculation was the
proposition of a higher order Laplacian model for both two and three (Equation
16) dimensional simulations [35, 54].

${\nabla{}}^{2}{\varphi{}}_{i}=\frac{1}{n_{0}}\sum_{i\not=j}\left(\frac{2{\varphi{}}_{ij}r_{e}}{r_{ij}}\right)$
(16)

where $\varphi$ is a generic physical quantity. This new derivation was named
by the authors as MPS with a Higher order Laplacian of pressure (MPS-HL).

####  3.3.3 Numerical stability

Khayyer and Gotoh [36] came up with a PPE’s source term with error-
compensating parts to enhance even further pressure and velocity field
calculations. The compensating parts should be measures for instantaneous and
accumulative violations of fluid incompressibility. Equation 17 shows the new
terms, and Equation 18 shows the complete modified pressure calculation
equation.

$\displaystyle
ECS=\left|{}\left(\frac{n^{k}-n_{0}}{n_{0}}\right)\right|{}\left[\frac{1}{n_{0}}{\left(\frac{Dn}{Dt}\right)}_{i}^{k}\right]+\left|{}\left(\frac{\Delta{}t}{n_{0}}{\left(\frac{Dn}{Dt}\right)}_{i}^{k}\right)\right|{}\left[\frac{1}{\Delta{}t}\frac{{n^{k}-n}_{0}}{n_{0}}\right]$
(17)

$\nabla{}{{}^{2}}p_{i}^{k+1}=\frac{\rho{}}{n_{0}\Delta{}t}{{\left(\frac{Dn}{Dt}\right)}_{i}}^{*}+ECS$
(18)

According to Gotoh [55], this method ensures satisfactory accuracy and stable
computation, more specifically, under the absence of tensile stress. Ref. [55]
shows a comparison between it and the standard MPS.

##  4 Implementation

This work implements the MPS through the C/C++ programming languages. OpenMP
[25] and CUDA [26] were used to take advantage of the many cores in the CPU
and GPU, hence accelerating the program execution. Lastly, there is a
graphical user interface (GUI) to aid in the system usage, developed through
the Windows Forms graphical class library from the Microsoft .NET framework
[56].

###  4.1 Neighbouring search algorithms

This study adopts the ”cell-linked list” strategy for the neighbourhood search
[43]. It relies on a background Cartesian grid that divides the whole domain
into cells. They have sides equal to $r_{e}$, which is the influence radius of
a particle. In every iteration, the particles of the simulation are allocated
in a specific cell, depending on their position. Thereby, when searching for a
particle’s neighbours, it is only necessary to look in the cell of the
particle itself and adjacent cells. Thus, there is a list of particles that
remain constant for that entire step, which offers a considerable decrease in
the complexity of the neighbouring search function to $O(n\ log\ n)$.

###  4.2 OpenMP

The machine used in the present study holds an Intel® Core™ i7-6820HK CPU @
2.70 GHz [57] with 32 GB of installed RAM, a 64-bit operating system (x64),
with 4 cores. 1 shows an example of the particle number density function main
loop parallelized just by adding one line of code, highlighted in red (line
1).

Listing 1: C code of the particle number density calculation benefiting from
OpenMP

⬇

1 #pragma omp parallel for schedule (guided)

for (int i = 1; i <= num_of_particles; i++) {

3 double sum = 0.0;

for (int l = 2; l <= neighb[i][1]; l++) {

5 int j = neighb[i][l];

(...) /*sum of kernel values for particle i and each of its neighbors j*/

7 }

n[i] = sum;

9}

return(n);

###  4.3 CUDA

The GPU code was developed based on the fully sequential and the OpenMP
versions previously presented using CUDA C/C++. The utilized machine has an
NVIDIA GeForce GTX 1080 Ti, which contains a total of $3,584$ CUDA cores and
11 GB of video memory. The number of threads per block is set to 256. As for
the implementations, while the parallelization process of some functions was
straightforward, some others needed adaptations for a parallelized version.

A tricky parallelization step in the fully incompressible model of the MPS is
the PPE solution. For this work, the standard ICCG solver used by Koshizuka
and his colleagues [45] was parallelized in both OpenMP and CUDA versions
using raw array formats for the linear system elements, such as the
coefficient matrix and right side source vector. Also, to decrease memory
usage by data structures, the total number of possible neighbours
($N_{neigh}$) of a particle was limited to $300$, which is the value used in
[8]. Doing so allowed assembling a coefficient matrix with $N\times N_{neigh}$
elements with $N$ being the total number of particles, rather than an $N\times
N$ matrix. 2 shows the PPE assembly in a CUDA kernel.

Listing 2: PPE assemble GPU code

⬇

1 unsigned int i = offset + (blockDim.x * blockIdx.x + threadIdx.x);

double val = 0.0;

3if (bcon[i] == 0) { // if fluid or boundary p.

poiss[i*num_of_neighbors + 1] = 0.0;

5 for (int l = 2; l <= neigh[i*num_of_neighbors + 1]; l++) {

int j = neigh[i*num_of_neighbors + l];

7 if (bcon[j] == -1) poiss[i*num_of_neighbors + l] = 0.0; //if dummy p.

else {

9 /* Calculating distance ‘dist’ between particles ... */

(...)

11 /*Calculating value ‘val’ for the matrix assemble ...*/

poiss[i*num_of_neighbors + l] = -val;

13 poiss[i*num_of_neighbors + 1] += val;

}

15 }

poiss[i*num_of_neighbors + 1] += (1.00e-7) / time_step / time_step;

17}

Another significant step of the method, performance-wise, is the neighbourhood
search function. This step occupies a large portion of the simulation running
time in either pressure calculation models, making its parallelization
essential to the practicability of massive simulations. In this work, the
implementation of the parallelized neighbourhood search is also based on the
cell linked list approach described in subsection 4.1, but, its
parallelization process is based on the one from Ref. [58], which presents a
high scalability capacity.

###  4.4 Graphical user interface

The user interaction may happen through a GUI. Figure 1 shows its layout and
design. In it, (1) is a combo box from which the user can choose whether the
simulation will have two or three dimensions; in (2) the user can choose how
the code will run: sequentially, in parallel in CPU (through OpenMP) or
parallel in GPU (through CUDA); in (3) the user will select a previously built
simulation scenario; in (4) two approaches of pressure calculation can be
chosen: weakly compressible or fully incompressible. By selecting the latter,
(5) and (6) will be available, which are options that, if checked, enable
Khayyer and Gotoh subsection 3.3 models of pressure gradient calculations to
better conserve the linear and angular momentum of the fluid flow. The
checking of (7) employs the SPS-LES turbulence model; (8) is only available if
the test scenario chosen allows multiphase interaction. If checked, it enables
viscoplastic properties in the second fluid in the simulation. (9) and (11)
are the density values [$Kg/m^{3}$] of the fluids in the simulation, and (10)
and (12) are their kinematic viscosity [$m^{2}/s$]; (13) sets the time
duration between two steps of the simulation and (14) sets how long it will
last, both in seconds. If checked, (15) creates a folder in the executable
file directory containing all the particles’ information in each time step. In
contrast, (16) creates a folder with the fluid particles’ information in each
time step. (17) starts generating the simulation, and (18) allows the user to
switch between languages.

Figure 1:  User interface main screen

##  5 Evaluation & discussion

The aim here is to validate the models and properties of the developed method,
such as incompressibility, numerical precision, turbulent flow, and multiphase
interaction. Another goal in this section is to assess the performance gain
and frame rates of the parallelized implementations of the weakly compressible
and fully incompressible versions of the MPS and discuss them.

###  5.1 Water drop

The examination of the evolution of an elliptic water drop is a common test
case to validate incompressibility models of a particle-based fluid simulation
method [47] [59] [43]. The usual test consists of a two-dimensional water
drop, beginning the simulation in the shape of a circle, with a predefined
velocity field of $(-100x,100y)~{}m/s$ so that its format evolves into an
elliptical shape over time. Figure 2 shows a sketch of the test’s geometry.
The water drop has a circle radius of $9\times 10^{-2}~{}m$ and an average
particle distance of $5\times 10^{-2}$ $m$, implying a total of $1257$
particles. The influence radius is $r_{e}=3l_{0}$, where $l_{0}$ is the
average particle distance, and the adopted time was $10^{-5}~{}s$. The used
configuration for this test was the default one, as presented in Figure 1,
with a fluid density equal to $10^{3}~{}Kg/m^{3}$ and fluid viscosity
$10^{-6}~{}m^{2}/s$, only differing in the field (3) since it is the Water
drop test.

Figure 2:  Initial state for the simulation in the water drop test (scale in
meters)

As in [43], Figure 3 shows three time instants of the simulation. The kernel
function used here is in Equation 4.

(a)  $t=0.0025~{}s$

(b) $t=0.0050~{}s$

(c) $t=0.0075~{}s$

Figure 3: Velocity magnitude profiles (in $m/s$) of the water drop in three
different time instants

As noted by Monaghan [47], the condition in this test to measure
incompressibility is that $ab$ is constant throughout the simulation, where
$a$ is the semi-minor axis, and $b$ is the semi-major axis of the ellipse.
This simulation ends when the size of $b$ becomes twice the value it was
initially, and, in that instant, the $ab$ value is compared to its initial
value. Errors are within less than $0.3~{}\%$. Table 1 compares this result to
[47] which keep errors in less than $2~{}\%$ and [43] which keep errors in
less than $0.2~{}\%$. This comparison shows the achieved result matches those
from previous works.

Table 1: Comparison of water drop dimensions throughout the simulation | $t=0.0~{}s$ | $t=0.01~{}s$ |  | Difference in %
---|---|---|---|---
$b$ (meters) | $0.875$ | $1.75$ |  | $100$
$ab$ (meters) | $0.766$ | $0.764$ |  | $0.26$

###  5.2 Dam break

The collapse of a water column has been widely used in the literature to
validate the numerical precision of various fluid simulation techniques.

As the test performed by [8], in this work, the water column height is two
times bigger than the water column length $L$. The floor in the model employed
here is also four times the length of the water column. The size of the water
column varies depending on how many particles the simulation has. The average
particle distance is $10^{-2}~{}m$, and the time step of the simulation is
$5\times 10^{-3}~{}s$. The total number of particles is $1122$.

The authors of the original MPS, referred here as standard MPS, put it to test
by comparing it to a volume of fluid (VOF) approach [60] and experimental data
from three different experiments [61, 62]. This test adopts the standard dam
break model [8] and compares the simulation results obtained with the
experimental data and the other simulation results. This comparison is
possible by examining the water leading-edge position over time, since the dam
burst until it hits the right wall. The leading edge is the front of the
collapsing water column running on the floor (bottom wall).

Figure 4 compares directly experimental data, other methods results, with the
proposed technique. The configuration used here is the default as well, as
depicted in Figure 1. The leading edge position is dimensionless, $Z/L$, where
$L$ is the water column’s initial length, which for this test is equal to
$0.18~{}m$. If the water column size changes, the value of $0.18~{}m$ also
changes, respecting the average particle distance $l_{0}$. The time axis in
the chart is dimensionless as well, $t\sqrt{2g/L}$, where $t$ is the time in
seconds and $g$ the gravitational acceleration, which is equal to
$9.8~{}m/s^{2}$.

Figure 4:  Direct comparison of the evolution of leading edge position between
the developed technique (in red), and the experimental results and other
methods. Background image extracted from [8]

Figure 4 shows that the developed method gets a lot closer to the experimental
data than the others. Interestingly, it almost overlaps the experiment
represented by the small empty squares from [61]. It is important to note
that, differently from a technique, the experiments do not become outdated
since the setup did not change over the years. Therefore, this corroborates
the assumption of a higher precision of the implemented method in this study.

###  5.3 R-T instability

The Rayleigh-Taylor instability test’s goal is to show the multiphase model
ability to handle the density stratification and to evolve a linear
perturbation into nonlinear hydrodynamic turbulence [63].

The test consists of placing the same amount of two immiscible fluids with
different densities and the same kinematic viscosity, one on the top of the
other. The heavier one will stay on top of the lighter, only influenced by the
gravitational force $g$ inside a two-dimensional rectangular box. When the
simulation starts, the denser fluid will tend to go downwards, pushing the
lighter upwards. The interface between the fluids will become unstable, and
the format of the pattern generated at the interface will say whether the
hydrodynamic turbulence was formed or not. The lighter fluid forms a bubble in
the shape of a mushroom cap that breaks eventually. The Atwood number,
calculated as shown in Equation 19, characterizes the problem.

$A_{t}=\frac{\rho_{h}-\rho_{l}}{\rho_{h}+\rho_{l}}$ (19)

$\rho_{h}$ and $\rho_{l}$ refer to the heavier and lighter fluid densities,
respectively. For the test used in this study, $A_{t}=1/3$ and the kinematic
viscosity $\nu=0.010$ for both fluids, following [44]. The test has $3066$
particles. The system configuration used in this scenario differs from the
default options shown in Figure 1 with some differences. Here, field (3) is
R-T instability, (4) is Weakly compressible, and fields (9) and (11) respect
the Atwood number calculation. 5(a) shows some steps of the generated
simulation. Ref. [44] shows a similar test displaying the same phenomenon. For
a better context on how other particle-based approaches behave, their results
are shown next to the generated simulation in 5(b) . The main difference
between them is that, in the scenario from Ref. [44], there is already an
initial perturbation to the fluids’ interface. Another difference is that in
the test performed here, the domain top is open. These differences influence
the test, causing the bubble (mushroom cap) to be upside-down.

(a) Rayleigh-Taylor instability generated by a multiphase interaction
developed in this work

(b) Rayleigh-Taylor instability by a basic multiphase MPS. Extracted from [44]

Figure 5: Rayleigh-Taylor instability in almost similar scenarios

The proposed technique generally shows good agreement regarding the Rayleigh-
Taylor instability when compared to another particle-based solution.

###  5.4 Oil spill

The Oil spill test makes possible a qualitative analysis and a unique
comparison to validate the multiphase model. [40] provides a replicable test
of a continuous oil spill due to a damaged tank.

The average particle distance is $4\times 10^{-3}~{}m$, the water and oil’s
kinematic viscosity $\nu_{water}$ and $\nu_{oil}$ are, respectively,
$10^{-6}~{}m^{2}/s$ and $5\times 10^{-5}~{}m^{2}/s$ and configuration here,
relative to Figure 1, differs in fields (3), which is Oil Spill, (4), Weakly
compressible, (11), which is $897.0Kg/m^{3}$ and (12), which is $5\times
10^{-5}~{}m^{2}/s$. This test has a total of $152,779$ particles.

This leakage simulation is qualitatively compared to experimental results and
to [40], which is a fully incompressible multiphase MPS. Figure 6 shows five
time instants of the experiments and the simulations. In each time instants,
there is: an image of the experiment; the simulation by [40]; and the
generated simulation, in that order.

(a) $t=0.2~{}s$

(b) $t=0.6~{}s$

(c) $t=3~{}s$

(d) $t=7~{}s$

Figure 6: Comparison between experiment, FS-MMPS [40] and the developed system
at different moments of the oil spill

The authors in [40] claim that FS-MMPS is quite precise and accurate. On the
other hand, it is possible to observe that the use of a weakly compressible
approach of the combined method developed here can better predict certain
situations, such as the more straight and even flow of the oil, in yellow,
over the denser fluid, mainly in 6(c) and 6(d). In a general manner, the oil
profile and the waves generated by it are in reasonable agreement with the
generated by the experiment, with fewer computations than the FS-MMPS, since,
in this case, a WC model is used, which significantly diminishes the pressure
calculation complexity.

###  5.5 Pressure field

It is essential to show how the previously presented numerical improvements
impact the pressure field of the fluid. This test makes use of the traditional
dam break scenario described in subsection 5.2. Figure 7 shows three time
instants of the same simulation. The left side refers to the simulation using
the standard MPS, and the right side, the simulation with the numerical
improvements described in subsection 3.3.

Figure 7:  Dam break simulation. Left side: standard MPS; Right side: modified
MPS with numerical improvements described in subsection 3.3. Time instants
referring 0.25 s, 0.5 s and 0.875 s into the simulation.

It is apparent in Figure 7 that the numerical improvements significantly
diminish the spurious pressure oscillations and can provide a more stable
pressure field during the simulation, even after the fluid suffers large
deformations.

###  5.6 Computational performance analysis

It is noticeable that switching between the incompressibility models - weakly
compressible (WC) or fully incompressible (FI) - had a material impact on
computational performance, both in memory usage and in runtime performance.
This discrepancy happens because the FI approach requires the solution of a
linear system of order $N\times N_{neigh}$, to ensure incompressibility and,
thus, obtain more precise values for the pressures. Storing it in memory and
solving it are both costly steps for the algorithm. The chosen
incompressibility model influences the most on the runtime, regardless of the
selected combination of enhancements.

To evaluate the computational performance, this section details the memory
usage and the runtime profile of two main approaches (WC and FI) combined with
three different types of execution: completely sequential; parallelized
through OpenMP; parallelized through CUDA. The entire runtime for each test
size was measured to calculate the speedups. Finally, it is possible to
calculate each version’s frame rate by running a set of tests in which they
vary in the number of particles.

####  5.6.1 Memory usage

The Performance and Diagnostics tool of Visual Studio 2019 [64] provides a
detailed report of the sequential and OpenMP executions to evaluate CPU memory
usage in all versions. The NVIDIA Visual profiler [65] provided the majority
of information regarding the GPU, such as the whole simulation runtime and
time spent in each CUDA kernel during the execution. As for the GPU memory
usage, GPU-Z [66] provided the minimum and maximum occupancy throughout the
simulation runtime.

Table 2 shows the memory usage of the WC approach of the MPS for a standard 2D
dam break test, like the one presented in subsection 5.2. For this case, there
is a total of $1,054,122$ particles.

Table 2: Occupancy in the CPU and GPU memories of different executions of the WC-MPS for 1 million particles | Sequential | OpenMP | CUDA
---|---|---|---
CPU Mem. usage | $4.7~{}\%$ | $4.7~{}\%$ | $4.4~{}\%$
GPU Mem. usage | 0 % | 0 % | $16.4~{}\%$

Table 3 shows the memory usage relative to the FI approach for the same test
described above and in subsection 5.2. One million particles can be considered
a high amount of particles for a FI approach running on GPU since a linear
system has to be fully loaded. This number of particles is an improvement
compared to previous works such as [67, 42].

Table 3: Occupancy in the CPU and GPU memories of different executions of the FI-MPS for 1 million particles | Sequential | OpenMP | CUDA
---|---|---|---
CPU Mem. usage | $20.3~{}\%$ | $20.3~{}\%$ | $20.0~{}\%$
GPU Mem. usage | 0 % | 0 % | $60.0~{}\%$

Note that, in both approaches, the CUDA versions require not only the GPU
memory but also CPU memory, which is also smaller compared to OpenMP and the
fully sequential executions. This difference may happen since some information
must exist in host memory (CPU) and in device memory (GPU) to transfer
particle input information from the CPU to the GPU.

####  5.6.2 Speedup

Due to the scalability of the GPU-accelerated neighbourhood search strategy,
the system was able to reach promising speedup values. Figure 8 shows that
this is especially true in more massive simulations. The OpenMP version did
not achieve considerable speedups since its neighbourhood calculation is not
optimally parallelized.

Figure 8:  Speedup caused by GPU parallelization versus total number of
particles in the simulation

Figure 8 shows that a case with $1,563,658$ particles can achieved a massive
speedup of $43.5$ times. The speedup curve shows that simulations with more
particles could yield even higher speedup rates. On the other hand, for a case
with $1,022$ particles, the speedup rate obtained was $0.8$ times, in other
words, running this scenario in the GPU is slower than running it in
sequentially in the CPU. That shows it is not efficient to employ parallel
programming for simulations with a small number of particles. One of the
causes for this is that the data transfer between the host and device reduces
the speedup significantly.

####  5.6.3 System limits

Data related to performance and memory limitations are presented here, such as
frame rate achieved for different numbers of particles and the maximum number
of particles in a simulation.

Figure 9 shows how many frames per second (FPS) the GPU-accelerated
application achieves, given the total number of particles in the simulation.
The blue frame rate curve represents a regular execution where the
neighbourhood is updated every step, prioritizing a more accurate simulation.
The red frame rate curve represents an execution where the list of neighbours
is only updated every four steps during the simulation, aiming for high
computational performance.

Figure 9:  Frame rate of the simulation given its total number of particles.

Hence, in a case that computational performance is the priority, Figure 9
shows that it is possible to simulate a scenario with approximately $200,000$
particles at one frame per second. It also displays a case with $46,768$
particles run at approximately 4 fps. Both of them are already considered
interactive frame rates [68, 69].

Regarding the simulation size limit, it is noteworthy that the available
memory plays a fundamental role. The GPU used in this work provides 11 GB of
video memory and allows for a simulation with $1,563,658$ before the
application runs out of memory. The largest simulation size tested for the
fully sequential version when adopting the fully incompressible approach had a
total of $3,677,442$ particles. The CPU process of the latter used 25 of the
32 GB available of RAM. When adopting the weakly compressible model, the total
number of particles could reach $16,657,962$ by using 26 GB of the available
RAM. The last two simulation loads show how expensive it is in terms of
memory, just loading the PPE into the RAM, not to mention the computational
load required to solve it.

##  6 Conclusions

As previously discussed, the study of fluid flow simulation is of great
importance in mitigating the consequences of environmental disasters and
accidents. It has applications in a wide range of engineering problems,
computer graphics, and virtual and augmented reality software. Meshless
methods like the MPS are a great alternative to deal with large deformations
and free-surface flow, situations where the traditional mesh-based approaches
usually perform inefficiently.

Throughout the development of his work, it became clear that the community is
continually improving the MPS, both in numerical performance and computation
efficiency, despite suffering from a few pressure instability problems. The
literature shows its usability in a large number of scenarios. The
considerable amount of referenced works also shows the complexity of this
task, the importance of the method, and the great potential to simulate,
increasingly more realistically, fluid flows.

MPS optimization is moderately complex since it is used to replicate real
phenomena more reliably, when taking into consideration other meshless
approaches. Some works enhance its computational efficiency with acceleration
structures without losing the precision it offers [23] [22] [67] [42]. In
contrast, other works prefer, despite the precision loss, replace the
performance and memory bottleneck of the method, the solution of a PPE, with
an equation of state to solve the pressures [43] [39] and then, parallelize it
through GPU [24] [70]. Another factor to consider is hardware development
since manufacturers are continuously building more powerful GPUs, which
directly influence computational performance.

The literature also shows that works usually accelerate the standard version
of the MPS [8] or the standard weakly compressible version [43][39], with few
numerical improvements to the calculation. This work provides a wide variety
of models, improvements, and approaches to the MPS technique, which are
entirely parallelized, through OpenMP and CUDA, integrated into a single
solution. This integration enables a fine-tuning of the system, allowing
setups concerned with high precision, fully incompressible fluid flow, and
GPU-accelerated multiphase WC fluid simulations. This solution is open-source
under the GPLv3 license.

Regarding the numerical improvements, the techniques proposed by [55] and [44]
were combined and extended, which implicated in enhancements in simulation
coherence, as presented in the Dam break, Oil spill, and Water drop tests. The
developed method shows compatibility with recent works in the Oil Spill test,
by qualitatively comparing it with [40], which shows the numerical advances
achieved. The GPU-acceleration provides speedups ranging from 3 to 43 times,
depending on the total number of particles in the simulation. That permits a
simulation with approximately $200,000$ particles to run at one frame per
second, and a test case with $46,768$ particles to achieve nearly 4 fps, which
both are already considered interactive rates [68, 69].

###  6.1 Future Work

There are exciting possibilities for future developments of this work.
Inevitably, refinements in every part of the code will lead to a more
optimized version, which is the path to a notable real-time simulation.

A possibility that would enhance even further computational performance is to
improve the GPU implementation of the FI-MPS to be able to run locally on
multiple GPU or in a GPU cluster, which could raise the speedup values to new
levels. [24] shows the performance of WC-MPS aided by a GPU cluster; fully
incompressible versions of the MPS could benefit significantly from such
structures. [71] presents a domain decomposition strategy for a parallelized
MPS for running in a cluster of computers. That enables massive simulations
since the simulation domain can be loaded and solved separately in the PC’s
memories and then integrated back.

Another possibility is to improve the OpenMP parallelization of the
neighbourhood search function. Since that function is based on its sequential
version, it could not reach decent speedup values. However, with an optimal
CPU parallelization of that step, the OpenMP acceleration could reach new
levels.

Still on the acceleration task, since solving a PPE is costly, an alternative
to the code parallelization could be the usage of neural networks. They would
learn the usual results of commonly yielded linear systems by fully
incompressible fluid simulations. The work of [72] shows promising fluid
simulations using this approach.

The addition of models that allow the interaction between fluids and solids,
such as floating bodies, deformable bodies, and viscoelastic fluid, would
drastically increase the number of application possibilities.

## Acknowledgements

The authors would like to thank Mr. Deep Tavker for authorizing the use of
part of his work to help develop the neighbourhood search strategy adopted
here.

## References

  * [1] PW Cleary, M Prakash, and J Ha.  Novel applications of smoothed particle hydrodynamics (SPH) in metal forming.  Journal of materials processing technology, 177(1):41–48, 2006\.
  * [2] Matthias Müller, David Charypar, and Markus Gross.  Particle-based fluid simulation for interactive applications.  In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’03, pages 154–159, Aire-la-Ville, Switzerland, Switzerland, 2003. Eurographics Association.
  * [3] Stefan Daenzer, Kevin Montgomery, Rüdiger Dillmann, and Roland Unterhinninghofen.  Real-time smoke and bleeding simulation in virtual surgery.  In MMVR, pages 94–99, 2007.
  * [4] Robert Bridson.  Fluid simulation for computer graphics.  CRC press, 2015.
  * [5] Ted Belytschko, Yury Krongauz, Daniel Organ, Mark Fleming, and Petr Krysl.  Meshless methods: an overview and recent developments.  Computer methods in applied mechanics and engineering, 139(1):3–47, 1996.
  * [6] Andrew A Johnson and Tayfun E Tezduyar.  Advanced mesh generation and update methods for 3d flow simulations.  Computational Mechanics, 23(2):130–143, 1999.
  * [7] Pascal-Jean Frey and Frédéric Alauzet.  Anisotropic mesh adaptation for CFD computations.  Computer methods in applied mechanics and engineering, 194(48):5068–5082, 2005.
  * [8] S Koshizuka and Y Oka.  Moving-particle semi-implicit method for fragmentation of incompressible fluid.  Nuclear science and engineering, 123(3):421–434, 1996.
  * [9] Robert A Gingold and Joseph J Monaghan.  Smoothed particle hydrodynamics: theory and application to non-spherical stars.  Monthly notices of the royal astronomical society, 181(3):375–389, 1977.
  * [10] Leon B Lucy.  A numerical approach to the testing of the fission hypothesis.  The astronomical journal, 82:1013–1024, 1977.
  * [11] Hitoshi Gotoh and Abbas Khayyer.  Current achievements and future perspectives for projection-based particle methods with applications in ocean engineering.  Journal of Ocean Engineering and Marine Energy, 2(3):251–278, Aug 2016.
  * [12] Hitoshi Gotoh and Tetsuo Sakai.  Lagrangian simulation of breaking waves using particle method.  Coastal Engineering Journal, 41(03n04):303–326, 1999.
  * [13] Abbas Khayyer.  Improved particle methods by refined models for free-surface fluid flows.  PhD thesis, Kyoto University, 2008.
  * [14] Hitoshi Gotoh, Hiroyuki Ikari, Tetsu Memita, and Tetsuo Sakai.  Lagrangian particle method for simulation of wave overtopping on a vertical seawall.  Coastal Engineering Journal, 47(02n03):157–181, 2005.
  * [15] Abbas Khayyer and Hitoshi Gotoh.  Modified moving particle semi-implicit methods for the prediction of 2d wave impact pressure.  Coastal Engineering, 56(4):419–440, 2009.
  * [16] Byung-Hyuk Lee, Jong-Chun Park, Moo-Hyun Kim, and Sung-Chul Hwang.  Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads.  Computer methods in applied mechanics and engineering, 200(9):1113–1125, 2011.
  * [17] Kazuya Shibata and Seiichi Koshizuka.  Numerical analysis of shipping water impact on a deck using a particle method.  Ocean Engineering, 34(3):585–593, 2007.
  * [18] Hitoshi Gotoh and Tetsuo Sakai.  Key issues in the particle method for computation of wave breaking.  Coastal Engineering, 53(2):171 – 179, 2006.  Coastal Hydrodynamics and Morphodynamics.
  * [19] Lei Fu and Yee-Chung Jin.  Investigation of non-deformable and deformable landslides using meshfree method.  Ocean Engineering, 109:192 – 206, 2015.
  * [20] Kazuya Shibata, Seiichi Koshizuka, Mikio Sakai, and Katsuji Tanizawa.  Lagrangian simulations of ship-wave interactions in rough seas.  Ocean Engineering, 42:13 – 25, 2012.
  * [21] Sung-Chul Hwang, Abbas Khayyer, Hitoshi Gotoh, and Jong-Chun Park.  Development of a fully lagrangian MPS-based coupled method for simulation of fluid-structure interaction problems.  Journal of Fluids and Structures, 50:497 – 511, 2014.
  * [22] XiaoSong Zhu, Liang Cheng, Lin Lu, and Bin Teng.  Implementation of the moving particle semi-implicit method on GPU.  SCIENCE CHINA Physics, Mechanics & Astronomy, 54(3):523–532, 2011\.
  * [23] Chiemi Hori, Hitoshi Gotoh, Hiroyuki Ikari, and Abbas Khayyer.  GPU-acceleration for moving particle semi-implicit method.  Computers & Fluids, 51(1):174–183, 2011.
  * [24] D Taniguchi, LM Sato, and LY Cheng.  Explicit moving particle simulation method on GPU clusters.  Blucher Mech. Eng. Proc. 1, 1:1155, 2014.
  * [25] L. Dagum and R. Menon.  Openmp: an industry standard api for shared-memory programming.  IEEE Computational Science and Engineering, 5(1):46–55, Jan 1998\.
  * [26] NVIDIA.  CUDA Zone — NVIDIA Developer, 2007.  Accessed: 2016-01-09.
  * [27] G. Chen, Y. Onishi, L. Zheng, and T. Sasaki.  Frontiers of Discontinuous Numerical Methods and Practical Simulations in Engineering and Disaster Prevention.  Taylor & Francis Group, London, 8 2013.
  * [28] M. Almeida, C. Brito, A. L. B. Vieira-e-Silva, V. Teichrieb, and J. M. Barbosa.  Meshless methods.  In G. Assi, H. Brinati, M. de Conti, and M. Szajnbok, editors, Applied Topics in Marine Hydrodynamics, chapter 8, pages 8.1–8.38. Escola Politécnica da Universidade de São Paulo (ISBN 978-85-86686-89-4), São Paulo, 2016.
  * [29] Fuji Technical Research Inc.  MPS-RYUJIN, 2013.  Accessed: 2018-04-23.
  * [30] Prometech Software.  Particleworks, 2014.  Accessed: 2018-04-23.
  * [31] Ahmad Shakibaeinia.  MPARS - Mesh-free Particle Simulator, 2012.  Accessed: 2018-04-09.
  * [32] Alejandro JC Crespo, José M Domínguez, Benedict D Rogers, Moncho Gómez-Gesteira, S Longshaw, R Canelas, Renato Vacondio, A Barreiro, and O García-Feal.  Dualsphysics: Open-source parallel CFD solver based on smoothed particle hydrodynamics (SPH).  Computer Physics Communications, 187:204–216, 2015.
  * [33] Alexis Hérault, Robert A. Dalrymple, Billy Edge, Giuseppe Bilotta, and Agnès Leroy.  GPUSPH.  Accessed: 2019-05-21.
  * [34] Masahiro Kondo, Kentaro Suto, Mikio Sakai, and Seiichi Koshizuka.  Incompressible free surface flow analysis using moving particle semi-implicit method.  Joint International Workshop: Nuclear Technology and Society – Needs for Next Generation, Berkeley, California, January 6-8, 2008, Berkeley Faculty Club, UC Berkeley Campus, 2008.
  * [35] Abbas Khayyer and Hitoshi Gotoh.  A higher order laplacian model for enhancement and stabilization of pressure calculation by the MPS method.  Applied Ocean Research, 32(1):124–131, 2010.
  * [36] Abbas Khayyer and Hitoshi Gotoh.  Enhancement of stability and accuracy of the moving particle semi-implicit method.  Journal of Computational Physics, 230(8):3093–3118, 2011.
  * [37] J.J. Monaghan.  SPH without a tensile instability.  Journal of Computational Physics, 159(2):290 – 311, 2000.
  * [38] Hitoshi Gotoh and Abbas Khayyer.  On the state-of-the-art of particle methods for coastal and ocean engineering.  Coastal Engineering Journal, 0(0):1–25, 2018.
  * [39] Ali Tayebi and Yee chung Jin.  Development of moving particle explicit (mpe) method for incompressible flows.  Computers & Fluids, 117:1 – 10, 2015.
  * [40] Guangtao Duan, Bin Chen, Ximin Zhang, and Yechun Wang.  A multiphase MPS solver for modeling multi-fluid interaction with free surface and its application in oil spill.  Computer Methods in Applied Mechanics and Engineering, 320:133 – 161, 2017.
  * [41] Eiji Harada, Hitoshi Gotoh, Hiroyuki Ikari, and Abbas Khayyer.  Numerical simulation for sediment transport using MPS-DEM coupling model.  Advances in Water Resources, 2017.
  * [42] André Luiz Buarque Vieira-e-Silva, Mozart W. S. Almeida, Caio Brito, and Veronica Teichrieb.  Improved MPS method and its variations for simulating incompressible fluids on GPU.  Journal on 3D Interactive Systems, 9(2), 2018.
  * [43] Ahmad Shakibaeinia and Yee‐Chung Jin.  A weakly compressible MPS method for modeling of open‐boundary free‐surface flow.  International Journal for Numerical Methods in Fluids, 63(10):1208–1232, 2010.
  * [44] Ahmad Shakibaeinia and Yee-Chung Jin.  MPS mesh-free particle method for multiphase flows.  Computer Methods in Applied Mechanics and Engineering, 229-232:13 – 26, 2012.
  * [45] Seiichi Koshizuka, Atsushi Nobe, and Yoshiaki Oka.  Numerical analysis of breaking waves using the moving particle semi-implicit method.  International Journal for Numerical Methods in Fluids, 26(7):751–769, 1998.
  * [46] G Batchelor.  K. 1967 an introduction to fluid dynamics, 1970.
  * [47] Joe J Monaghan.  Simulating free surface flows with SPH.  Journal of computational physics, 110(2):399–406, 1994.
  * [48] J J Monaghan.  Smoothed particle hydrodynamics.  Reports on Progress in Physics, 68(8):1703, 2005.
  * [49] Cx K Batchelor and GK Batchelor.  An introduction to fluid dynamics.  Cambridge university press, 2000.
  * [50] Songdong Shao and Hitoshi Gotoh.  Turbulence particle models for tracking free surfaces.  Journal of Hydraulic Research, 43(3):276–289, 2005.
  * [51] J. Smagorinsky.  General circulation experiments with the primitive equations.  Monthly Weather Review, 91(3):99–164, 1963.
  * [52] Robert S Rogallo and Parviz Moin.  Numerical simulation of turbulent flows.  Annual Review of Fluid Mechanics, 16(1):99–137, 1984.
  * [53] Abbas Khayyer and Hitoshi Gotoh.  Development of CMPS method for accurate water-surface tracking in breaking waves.  Coastal Engineering Journal, 50(02):179–207, 2008.
  * [54] Abbas Khayyer and Hitoshi Gotoh.  A 3d higher order laplacian model for enhancement and stabilization of pressure calculation in 3d MPS-based simulations.  Applied Ocean Research, 37:120–126, 2012.
  * [55] H Gotoh.  Advanced particle methods for accurate and stable computation of fluid flows.  Frontiers of Discontinuous Numerical Methods and Practical Simulations in Engineering and Disaster Prevention, page 113, 2013.
  * [56] David S. Platt.  Introducing Microsoft .Net, Second Edition.  Microsoft Press, Redmond, WA, USA, 2nd edition, 2002.
  * [57] Intel.  Intel processor i7 4790 specifications, 2013.  Accessed: 2018-04-22.
  * [58] Deep Tavker.  Parallel Neighbour Search Implementation.  https://github.com/deeptavker/Parallel-Neighbour-Search, 2018.
  * [59] J Bonet and T-SL Lok.  Variational and momentum preservation aspects of smooth particle hydrodynamic formulations.  Computer Methods in applied mechanics and engineering, 180(1):97–115, 1999.
  * [60] C.W Hirt and B.D Nichols.  Volume of fluid (vof) method for the dynamics of free boundaries.  Journal of Computational Physics, 39(1):201 – 225, 1981.
  * [61] Seichii Koshizuka.  A particle method for incompressible viscous flow with fluid fragmentation.  Comput. Fluid Dynamics J., 4:29–46, 1995.
  * [62] J. C. Martin, W. J. Moyce, J. C. Martin, W. J. Moyce, W. G. Penney, F. R. S., A. T. Price, and C. K. Thornhill.  Part iv. an experimental study of the collapse of liquid columns on a rigid horizontal plane.  Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 244(882):312–324, 1952.
  * [63] David L. Youngs.  Numerical simulation of turbulent mixing by rayleigh-taylor instability.  Physica D: Nonlinear Phenomena, 12(1):32 – 44, 1984.
  * [64] Microsoft.  Visual Studio Enterprise 2015, 2018.  Accessed: 2018-05-01.
  * [65] NVIDIA.  NVIDIA Visual Profiler, 2018.  Accessed: 2018-05-01.
  * [66] TechPowerUp.  GPU-Z Video card GPU Information Utility, 2018.  Accessed: 2018-05-01.
  * [67] A. L. B. Vieira-e-Silva, M. W. S. Almeida, C. Brito, and V. Teichrieb.  Improved meshless method for simulating incompressible fluids on GPU.  In 2017 19th Symposium on Virtual and Augmented Reality (SVR), pages 297–308, Nov 2017.
  * [68] Ingo Wald, Andreas Dietrich, and Philipp Slusallek.  An interactive out-of-core rendering framework for visualizing massively complex models.  In SIGGRAPH ’05, 2004.
  * [69] D. Overby, Z. Melek, and J. Keyser.  Interactive physically-based cloud simulation.  In 10th Pacific Conference on Computer Graphics and Applications, 2002. Proceedings., pages 469–470, 2002.
  * [70] Davi T. Fernandes.  Implementação de framework computacional de paralelização híbrida do Moving Particle Semi-implicit Method para modelagem de fluidos incompressíveis.  PhD thesis, Universidade de São Paulo, 2013.
  * [71] Davi Teodoro Fernandes, Liang-Yee Cheng, Eric Henrique Favero, and Kazuo Nishimoto.  A domain decomposition strategy for hybrid parallelization of moving particle semi-implicit (MPS) method for computer cluster.  Cluster Computing, 18(4):1363–1377, Dec 2015.
  * [72] Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin.  Accelerating eulerian fluid simulation with convolutional networks.  arXiv preprint arXiv:1607.03597, 2016.


# Finding the maximum-a-posteriori behaviour of agents in an agent-based model

Daniel Tang
Leeds Institute for Data Analytics
University of Leeds
Leeds, UK
<EMAIL_ADDRESS>
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. 757455)

###### Abstract

In this paper we consider the problem of finding the most probable set of
events that could have led to a set of partial, noisy observations of some
dynamical system. In particular, we consider the case of a dynamical system
that is a (possibly stochastic) time-stepping agent-based model with a
discrete state space, the (possibly noisy) observations are the number of
agents that have some given property and the events we’re interested in are
the decisions made by the agents (their “expressed behaviours”) as the model
evolves.

We show that this problem can be reduced to an integer linear programming
problem which can subsequently be solved numerically using a standard branch-
and-cut algorithm. We describe two implementations, an “offline” algorithm
that finds the maximum-a-posteriori expressed behaviours given a set of
observations over a finite time window, and an “online” algorithm that
incrementally builds a feasible set of behaviours from a stream of
observations that may have no natural beginning or end.

We demonstrate both algorithms on a spatial predator-prey model on a 32x32
grid with an initial population of 100 agents.

_K_ eywords Data assimilation, Bayesian inference, Agent based model, Integer
linear programming, predator prey model

##  1 Introduction

The problem we’ll consider here is as follows: there is a dynamical system of
interest which we believe can be well modelled by an agent based model. We
have some prior beliefs about the boundary conditions of the system and have
made some observations over some period of time. However, the observations are
incomplete (in that they do not specify the full state of the model at any
time) and possibly noisy. The problem is to find the most probable sequence of
behavioural choices made by the agents (their “expressed behaviours”) given
the observations. The expressed behaviours consist of a feasible set of agent
actions or events that could have given rise to the observations.

Before presenting an algorithm to solve this problem, we’ll first define the
problem a bit more precisely.

###  1.1 What is an agent-based model?

We take the essential property of an agent-based model to be that the dynamics
of the model can be described entirely in terms of agent behaviours (if the
environment contains objects, these should also be modelled as agents). Each
agent may have an internal state, and the state of the whole model at any
instant can be completely specified by the states of all the agents.

Within this broad definition, there are a number of sub-categorizations

  * •

Modelled time can be either discrete or continuous. In a discrete-time model,
agent behaviour consists of timesteps from time $t$ to $t+1$ and time is an
integer. In a continuous-time model agent behaviour consists of discrete
events that can happen at any time expressed as a real number. In this paper
we will consider discrete-time models. We will also suppose that every agent
gets an opportunity to act in each timestep, and that each agent potentially
has access to the states of other agents at time $t$ when making the timestep
to $t+1$.

  * •

The internal state of an agent can be discrete, continuous or mixed depending
on whether the state can be most naturally expressed as a vector of integers,
real numbers or a mixture of both. In this paper we will consider agents with
discrete-state.

  * •

Agents can be deterministic or stochastic. A deterministic agent will always
exhibit the same behaviour if its own internal state and the states of other
agents in the model are fixed, whereas a stochastic agent will exhibit
different behaviours randomly drawn from a fixed distribution under the same
circumstances. In this paper we will consider stochastic agents.

  * •

If only one agent can occupy a given state at any instant, then the state of
the model can be described as a set of agent states and we call it a Fermi-
Dirac model, whereas if two or more agents can occupy the same state at the
same time then the model state must be described as a multisetBlizard . (1988)
of agent states and we call it a Bose-Einstein model. In this paper we will
consider Bose-Einstein models.

We define the behaviour of an agent as a computer program which describes a
timestep of the agent. In addition to the usual programming constructs, the
program may:

  * •

query the state of the acting agent at time $t$

  * •

query whether there exist any other agents in the model at time $t$ that
satisfy any given predicate on their state

  * •

create a new agent in any state

  * •

remove the acting agent from the model

  * •

modify the state of the acting agent

  * •

call a function that returns a random number.

A timestep of the whole model consists of the execution of the timestepping
program for each agent.

Such a program defines a conditional probability distribution over a set of
events Staton . (2016). Each event represents an agent performing an action
that expresses a specific behaviour, which we will write as
$(\Phi|\psi,R,\bar{R})$ where $\psi$ is the internal state of the agent before
the behaviour is expressed, $R$ is the set of other agents (states) that must
be present for the behaviour to be expressed, $\bar{R}$ is the set of agents
that must be absent for the behaviour to be expressed and $\Phi$ is the
multiset of agents that are the result of the action occurring (including the
final state of the acting agent). So, if we say an agent has a behaviour
$P(\Phi|\psi,R,\bar{R})$ then when it is in state $\psi$ in an environment
containing all agents in $R$ and no agents in $\bar{R}$, it performs an action
that results in $\Phi$ with probability $P(\Phi|\psi,R,\bar{R})$.

So, for example, if an agent in state $\psi$ at time $t$ has a behaviour
$P(\\{\psi^{\prime}\\}|\psi,\emptyset,\emptyset)=1.0$ it will, with certainty,
change its state to $\psi^{\prime}$ at time $t+1$ . An agent with behaviour
$P(\emptyset|\psi,\\{\phi\\},\emptyset)=0.5$ in state $\psi$ in the presence
of another agent in state $\phi$ will disappear with probability $0.5$.

Without loss of generality, we assume that all agents in the model have the
same behaviour. If necessary, multiple types of agent can be modelled by
making the agent’s type part of its state.

In any timestep of the model a multiset of events will occur, we will call
such a multiset a model event to distinguish it from the agent events that are
members of the model event. The consequence of a model event, $E$, is the sum
of the consequences of its members. We’ll write this as $\Psi(E)$

$\Psi(E)=\sum_{(\Phi|\psi,R,\bar{R})\in E}\Phi$ (1)

where a sum over multisets is defined in the obvious way. We define an event
trajectory, $\mathcal{E}$, as a list of model events
$\left<E_{1}...E_{n}\right>$, which can be thought of as the assertion that,
for each $i$, the events in $E_{i}$ occurred in timestep from $t=i-1$ to
$t=i$. Not all event trajectories are feasible as, for each agent event, the
presence and absence conditions of the event must be met.

Given an event trajectory, $\mathcal{E}=\left<E_{1}...E_{n}\right>$, and a
seed agent, $\psi_{\emptyset}$, let

$\Psi_{i}=\begin{cases}\\{\psi_{\emptyset}\\}&\text{if }i=0\\\
\Psi(E_{i})&\text{otherwise}\\\ \end{cases}$

and let $\mathcal{E}$ be feasible if and only if $\forall i\in 1...n:$

Each agent must perform exactly one action per timestep

$\sum_{(\Phi|\psi,R,\bar{R})\in E_{i}}\\{\psi\\}=\Psi_{i-1}$ (2)

…each event’s presence-requirements must be fulfilled

$\bigcup_{(\Phi|\psi,R,\bar{R})\in E_{i}}R\subseteq\Psi_{i-1}$

…and each event’s absence-requirements must be fulfilled

$\left(\bigcup_{(\Phi|\psi,R,\bar{R})\in
E_{i}}\bar{R}\right)\cap\Psi_{i-1}=\emptyset$

The probability of a feasible model event is given by

$P(E)=\prod_{(\Phi|\psi,R,\bar{R})\in E}P(\Phi|\psi,R,\bar{R})$

and the probability of an event trajectory is just the product of the
probabilities of all the events

$P(\left<E_{0}...E_{m}\right>)=\prod_{i}P(E_{i})$

###  1.2 What is an observation?

In this paper we will consider observations that consist of a count of the
number of agents for which some predicate is true. The observation can be
exact or noisy in that the real count can lie in some range. More formally
$\omega=\left<L,U,B\right>$ denotes the observation that the number of agents,
$N$, for which predicate $B$ is true is somewhere in the range $L\leq N\leq
U$. If $E$ is a model event, then we say that $E$ satisfies the observation,
which we write as $E\vdash\omega_{i}$, if and only if

$L\leq\left|\left\\{\psi:\psi\in\Psi(E)\wedge B(\psi)\right\\}\right|\leq U$

From this, we say that $E$ satisfies a set of observations, $\Omega$, if it
satisfies all of its members. We can extend this to model event trajectories
$\mathcal{E}=\left<E_{1}...E_{n}\right>$ and say that $\mathcal{E}$ satisfies
a list of sets of observations (which we’ll call an observation of a
trajectory), $\mathcal{O}=\left<\Omega_{1}...\Omega_{n}\right>$ if for each
$i$, $E_{i}$ satisfies $\Omega_{i}$.

This form of observation can deal with exact observations ($L=U$),
observations where there may be false-positives ($U$ is the observed count and
$L=0$), observations where there may be false negatives ($L$ is the observed
count and $U=\infty$) or observations that have some uncertainty range
($L=N-\Delta N$ and $U=N+\Delta N$).

###  1.3 What is a prior belief about boundary conditions?

The boundary conditions on the system consist of agents entering the system at
particular times. At its simplest, it is a set of agents that enter the system
at the first timestep and are created by the seed agent, $\psi_{\emptyset}$,
through a behaviour $P(\Phi|\psi_{\emptyset},\emptyset,\emptyset)=1.0$, where
$\psi_{\emptyset}\notin\Phi$. More generally, it is a set of probabilities
$P(\Gamma,t)$ which give the prior probability that a multiset $\Gamma$ is a
subset of the agents that enter the system at time $t$.

In this paper we assume that these prior beliefs can be split into factors
that involve single agent states. In this case the prior probabilities can be
expressed as

$P(\Gamma,t)=\prod_{\psi\in^{=n}\Gamma}P(\psi,n,t)$

where $\psi\in^{=n}\Gamma$ is true when $\psi$ appears in $\Gamma$ exactly $n$
times and $P(\psi,n,t)$ is the probability of at least $n$ agents in state
$\psi$ being injected at time $t$.

This can conveniently be represented in terms of the behaviour of a set of
hypothetical (though never explicitly represented computationally) zero energy
agents $\emptyset_{\psi}$ which, are created by the seed agent
$\psi_{\emptyset}$. The behaviour of a zero energy agent has the form,
$(\\{\emptyset_{\psi},\psi^{n}\\}|\emptyset_{\psi},\emptyset,\emptyset)$ which
results in $n$ agents in the state $\psi$ being injected into the system,
along with the continued existence of the injecting, zero energy agent. If we
allow the behaviour of these agents to be time dependent (either implicitly or
explicitly via an internal state), then any set of boundary conditions can be
represented within the existing formalism.

##  2 An algorithm to find the MAP event trajectory

Bringing all this together, we can now express our problem precisely:

Given an observation of a trajectory $\mathcal{O}$ and an agent behaviour
$P(\Phi|\psi,R,\bar{R})$, find

$\mathcal{E}=\arg\max_{\left<E_{0}...E_{m}\right>}\prod_{i}\prod_{\epsilon\in
E_{i},}P(\epsilon)$ (3)

subject to $\mathcal{E}$ being feasible and satisfying $\mathcal{O}$.

To turn this into a tractable algorithm, first represent an event trajectory
as a set of non-negative integers $c_{te}$ such that in timestep $t$, the
event $e$ occurs exactly $c_{te}$ times.

Now encode each event in an agent’s behaviour, $e=(\Phi|\psi,R,\bar{R})$, by
letting

$\Phi_{e\phi}=n:\phi\in^{n}\Phi\\\ $ (4)

$\psi_{e\phi}=\begin{cases}1&\text{if }\phi=\psi\\\ 0&\text{otherwise}\\\
\end{cases}$ (5)

$R_{e\phi}=\begin{cases}1&\text{if }\phi\in R\\\ 0&\text{otherwise}\\\
\end{cases}$ (6)

$\bar{R}_{e\phi}=\begin{cases}1&\text{if }\phi\in\bar{R}\\\
0&\text{otherwise}\\\ \end{cases}$ (7)

Now encode the observation predicates for each event’s consequences

$B_{e}=|\\{\psi:\psi\in\Phi\wedge B(\psi)\\}|$ (8)

and let the probability of each event be written

$P_{e}=P(\Phi|\psi,R,\bar{R})$

finally, let the consequences of a model event be encoded as

$\Psi_{t\phi}=\sum_{e}\Phi_{e\phi}c_{te}$ (9)

and define a set of boolean variables $0\leq b_{t\phi}\leq 1$ which satisfy

$0\leq Mb_{t\phi}-\Psi_{t\phi}$ (10) $0\leq\Psi_{t\phi}-b_{t\phi}$ (11)

for some multiplicity $M$ which can be any number larger than the maximum
number of agents in any one state at any one time. This, for example, could be
set to the expected maximum number of agents in the model at any one time.

We can now express the problem as:

Maximise

$P=\sum_{t}\sum_{e}\log(P_{e})c_{te}$ (12)

subject to the observations being satisfied

$\forall\Omega_{t}\in\mathcal{O},\left<L,U,B\right>\in\Omega_{t}:L\leq\sum_{e}B_{e}c_{te}\leq
U$ (13)

…each agent performing exactly one action

$\forall\phi:0\leq\Psi_{(t-1)\phi}-\sum_{e}\ \psi_{e\phi}c_{te}\leq 0$ (14)

…each event’s presence requirements being fulfilled

$\forall\phi:0\leq Mb_{(t-1)\phi}-\sum_{e}\ R_{e\phi}c_{te}$ (15)

…each event’s absence requirements being fulfilled

$\forall\phi:\sum_{e}\ \bar{R}_{e\phi}c_{te}+Mb_{(t-1)\phi}\leq M$ (16)

This can be seen to be in the form of a pure integer linear programming
problem, which in full generality is NP complete, but for which there exists
algorithms that, in practice, are able to find solutions to problems of sizes
large enough to be useful.

##  3 Online assimilation

The algorithm described above assumes that the observations span a finite time
window. If the observations span a very long window or originate from a stream
which has no specific start or end, we would usually split the stream up into
smaller windows and deal with the observations one window at a time. However,
if we commit to the MAP event trajectory of one window, and use the end state
of this as the start state of the next window, we are in danger of being led
down “garden paths”; i.e. the observations from one window lead us into
believing the most probable explanation, only for later observations to show
us that a much less probable and completely incompatible explanation is in-
fact the case. So, we can’t commit to believing an event trajectory for one
window because, in the light of the observations from the next window, the
trajectory may turn out to be impossible. If this happens the assimilation
must stop because there is no event trajectory that satisfies the new window’s
observations while being a feasible extension of the previous window’s
trajectory. To solve this problem, we present a second algorithm for use with
longer windows and streaming observations.

We begin by defining the idea of a partial event trajectory which is like a
normal event trajectory except that instead of each agent having to perform
exactly one action per timestep, they can perform at most one action per
timestep. So, in place of equation 2 we have

$\sum_{(\Phi|\psi,R,\bar{R})\in E_{i}}\\{\psi\\}\subseteq\Psi(E_{i-1})$ (17)

Intuitively, a partial event trajectory can be interpreted as saying that at
least these events happened in these timesteps, but other additional events
may have happened concomitantly.

In an partial trajectory, all observations up to the end of the current window
are satisfied but only a subset of the agents present at any particular time
will be moved forward in time by the next model event, so agents may be left
hanging at a time in the middle of the trajectory. This is good because it
reduces the likelihood that we will be led down a garden path by delaying any
commitment to the exhibited behaviour of the hanging agents until more
observations are available.

The assimilation proceeds by finding a partial trajectory for one window, then
taking the next window of observations and finding the MAP partial trajectory
given these new observations along with the partial trajectory we’ve already
committed to. In order to correctly move hanging agents forward in time, we
need to distinguish between the events that we’ve already committed to and the
events we’re optimising over. If we let $\hat{c}_{te}$ be the number of events
of type $e$ that we’ve committed to in timestep $t$, then equation 18 becomes

$\Psi_{t\phi}-\sum_{e}\Phi_{e\phi}c_{te}=\sum_{e}\Phi_{e\phi}\hat{c}_{te}$
(18)

the agency constraint, equation 14, becomes

$\forall\phi:\sum_{e}\ \psi_{e\phi}\hat{c}_{te}\leq\Psi_{(t-1)\phi}-\sum_{e}\
\psi_{e\phi}c_{te}$ (19)

and the absence constraint, equation 16, becomes

$\forall\phi:\sum_{e}\ \bar{R}_{e\phi}c_{te}+Mb_{(t-1)\phi}\leq M-\sum_{e}\
\bar{R}_{e\phi}\hat{c}_{te}$ (20)

Since observations are fully satisfied by a partial trajectory, the
observation constraints, equation 13, now only apply to timesteps for which
there are no committed events. The presence requirement constraints can remain
unchanged.

Even though partial trajectories reduce the probability of us being led down
the garden path, they do not eradicate it completely. So, if we end up with
observations that are incompatible with our commitments, we rollback the
commitments, one timestep at a time, from the most recent until a solution
becomes possible. In the worst case, we’ll have to roll-back all the way to
$t=0$, but in practice the probability of this will be very small for most
models.

One final issue with online assimilation is agents that leave the model. When
an agent leaves the model it will not be observed again, but the online
assimilation will just leave it in the hanging state on the off-chance that
it’s still around but has evaded observation. At some point we are likely to
want to commit to saying that the agent has left the model. To do this we can
estimate the probability of the agent evading observation for so long, and
when this estimate falls below a certain level, we commit to its departure.

If required, a partial event trajectory can be extended into a complete
trajectory at any point by requiring that equation 17 is satisfied with
equality and solving the resulting set of equations.

##  4 Experiments with a spatial predator-prey model

In order to demonstrate the above algorithms, we simulated a spatial predator-
prey model on a 32x32 grid with periodic boundary conditions. Agents can be
either predator or prey. At each timestep, agents can stay on the same grid-
square, move to an adjacent grid-square (i.e. up, down, left or right) or die.
Prey can also give birth to another prey on an adjacent grid-square. Prey
cannot move onto a square that contains a predator, but in the presence of a
predator on an adjacent square they may be eaten and the predator reproduce
into the newly vacant square. The probability of each behaviour is shown in
table 1. Perhaps surprisingly, even a model with such a simple set of
behaviours exhibits quite complex emergent properties that are difficult to
predict.

Aget type | Description | Probability
---|---|---
Prey | die | 0.03
| reproduce | 0.06
| move (up/down/left/right) / be eaten | 0.728
| stay put / be eaten | 0.182
Predator | die | 0.05
| move (up/down/left/right) | 0.76
| stay put | 0.19
Table 1: The rates of each behaviour in the predator-prey model

Forward simulations of the ABM were performed in order to create synthetic
observation data. The initial state of each simulation consisted of 40
predators and 60 prey placed at random positions on the grid with uniform
probability. The boundary conditions for the MAP analyses consisted of a fully
observed initial state. Observations were made at each timestep and each agent
was observed with a probability $\frac{2}{3}$ in each timestep.

###  4.1 Offline MAP experiment

Using the above boundary conditions and observations and using a window length
of seven timesteps, we set up an integer linear programming problem as
specified above and used Google’s OR-ToolsGoogle (2020) library to find the
MAP event trajectory. In total 16 simulations were analysed, the mean time
taken per simulation to compute the MAP using a single core of a 1.6GHz Intel
i5-8250U was 156s.

A comparison between typical real and MAP model states at the end of an
assimilation window is shown in figure 1. As a rough measure of the difference
between the MAP and the “real” state, Figure 2 shows the average shortest
$L_{1}$ distance between the unobserved agents in the MAP and the nearest
unobserved agent of the same type in the real state for each timestep in the
window. For comparison, we also show the average distance for a model that
places the agents randomly with uniform probability.

Figure 1: Snapshot comparing the real and MAP states at the end of a 7-step
assimilation window. Blue symbols indicate the spatial position of predators,
red symols indicate prey. Squares show the position of agents in the real
state, while triangles show agents in the MAP state. Figure 2: The average
shortest $L_{1}$ distance from an unobserved agent in a MAP state to an
unobserved agent of the same type in the corresponding real state, plotted
against number of timesteps from the start of the window. For comparison, we
also show the shortest distance for the trajectories produced by the online
algorithm, and the shortest distance for a randomly placed agent.

###  4.2 Online assimilation experiment

The online algorithm was applied to the same set of observations and boundary
conditions as the offline experiments. A window size of one timestep was used
and agents were not removed from the model if unobserved. Once the partial
event trajectory for a whole simulation was found, the trajectory was
completed by solving one more time for the hanging agents. The mean time per
simulation to compute the event trajectory, using the same computer as used
for the offline experiments, was 11.4s, representing a speed increase of just
under 14 times compared to the offline algorithm.

In order to assess the quality of the online algorithm compared to that of the
MAP, we calculated the ensemble average log ratios of the posterior
probabilities

$\overline{\ln\left(\frac{P(\mathcal{E}_{MAP}|\mathcal{O})}{P(\mathcal{E}_{real}|\mathcal{O})}\right)}=66.7$

$\overline{\ln\left(\frac{P(\mathcal{E}_{online}|\mathcal{O})}{P(\mathcal{E}_{real}|\mathcal{O})}\right)}=44.1$

where $\mathcal{O}$ are the observations, $\mathcal{E}_{MAP}$ is the MAP event
trajectory, $\mathcal{E}_{online}$ is the trajectory given by the online
algorithm and $\mathcal{E}_{real}$ is the trajectory used to create the
observations. The MAP log ratio can be thought of as an upper bound on
performance while a log ratio of around 0 would be the expected value if we
sampled from the posterior. The performance of the online algorithm can be
assessed in relation to these two reference points. As can be seen, the online
algorithm performs less well than the MAP algorithm, but the online
trajectories are, on average, substantially more probable than the real
trajectory that was used to create the observations in the first place.

Figure 2 shows the average shortest $L_{1}$ distances between agents in the
the online trajectory and those in the real trajectory, alongside those of the
MAP trajectory.

###  4.3 Code

The code for all experiments described here is available at
https://github.com/danftang/MaxAPosteriori

##  5 Conclusion

We have demonstrated that the maximum-a-posteriori trajectory of an agent-
based, spatial predator-prey model can be found given a set of partial
observations. We have also shown how to construct a feasible interpretation of
a stream of such observations of arbitrary length.

Since integer linear programming is NP-complete the technique described here
will not in full generality scale to very large models but will work best when
agents choose from one of a limited number of options at each timestep, when
interactions are sparse and when the observations impose strong constraints on
the state. However, expressing the problem in this form lays the foundations
for further research to apply approximations that may lead to more scalable
algorithms.

Finally, it is worth noting that while the MAP event trajectory is the most
probable set of events given the observations, it is also usually highly
unlikely to be the true cause of the observations, for the simple fact that it
is one among a very large number of other trajectories with non-trivial
posterior probability. For this reason it is not usually a good idea to infer
values of unobserved variables from the MAP. Rather, the MAP should be used as
one component of a larger analysis. For example, the MAP trajectory is often a
good choice for the start state of a Markov chain Monte Carlo sampling
algorithm, which will be described in our next paper.

## References

  * Blizard . (1988) blizard1988multisetBlizard, WD. .  1988\.  Multiset theory. Multiset theory.  Notre Dame Journal of formal logic30136–66.
  * Google (2020) googleortoolsGoogle.  2020\.  OR-Tools. Or-tools.  https://developers.google.com/optimization [Online; accessed 18-February-2020]
  * Staton . (2016) staton2016semanticsStaton, S., Wood, F., Yang, H., Heunen, C. Kammar, O.  2016\.  Semantics for probabilistic programming: higher-order functions, continuous distributions, and soft constraints Semantics for probabilistic programming: higher-order functions, continuous distributions, and soft constraints.  2016 31st Annual ACM/IEEE Symposium on Logic in Computer Science (LICS) 2016 31st annual acm/ieee symposium on logic in computer science (lics) ( 1–10).