WATER ANIMATION USING COUPLED SPH AND WAVE

EQUATIONS
by
Varun Ramakrishnan

A Thesis
Submitted to the Faculty of Purdue University
In Partial Fulfillment of the Requirements for the degree of

Master of Science

Department of Computer Graphics Technology

West Lafayette, Indiana
May 2023
THE PURDUE UNIVERSITY GRADUATE SCHOOL
STATEMENT OF COMMITTEE APPROVAL

Dr. Tim McGraw, Chair

Department of Computer Graphics Technology

Dr. Christos Mousas

Department of Computer Graphics Technology

Dr. Yingjie Chen

Department of Computer Graphics Technology

Approved by:
Dr. Vetria Byrd

2
 ACKNOWLEDGMENTS

Im extremely grateful to my advisor—Prof. Tim McGraw, for his invaluable patience,

feedback, knowledge, and expertise. I also could not have undertaken this journey without
the rest of my committee—Prof. Christos Mousas, and Prof. Yingjie Chen, who provided
great guidance, support, and feedback.

3
 TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

LIST OF SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1 Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 Significance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5 Hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.6 Project Deliverables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Delimitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.9 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 LITERATURE REVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1 Water Animation using SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Water Animation using Wave Equations . . . . . . . . . . . . . . . . . . . . 22
2.3 Other Methods of Water Animation . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Ships, Splashes, and Waves on a Vast Ocean . . . . . . . . . . . . . . 23
2.3.2 Water Surface Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.3 Real-Time Cartoon Water Animation . . . . . . . . . . . . . . . . . . 24
2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 METHODOLOGY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1 Development Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4
3.2 Research Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.2 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.3 Acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.4 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.5 Kernel Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.3.6 Neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Wave Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.5 Program Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.1 Data Management . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Compute Shader . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Work Groups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Shader Storage Buffer Objects . . . . . . . . . . . . . . . . . . . . . . 37
3.6 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.6.1 Uniform Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
Collision Query . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
SSBO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6.2 SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Foam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Force Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Density and Pressure Computation . . . . . . . . . . . . . . . . . . . 46
Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
SSBO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.6.3 Wave Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.6.4 Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Splashes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
Breaking Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5
 Boat Wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.7 Rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.7.1 Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7.2 Wave Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.8.1 Visuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.8.2 Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.8.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4 SUMMARY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

6
 LIST OF TABLES

3.1 Memory needed for different particle counts . . . . . . . . . . . . . . . . . . . . 35

3.2 Cell index, count, and start arrays for grid in figure 3.2 . . . . . . . . . . . . . . 42
3.3 content list for grid in figure 3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Memory occupied by contents array for different particle counts . . . . . . . . . 44
3.5 Memory occupied by count, and start for different grid counts . . . . . . . . . . 44
3.6 Visual comparison of animations . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.7 Uniform Grid compute times in this study . . . . . . . . . . . . . . . . . . . . . 57
3.8 SPH compute times in this study . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.9 Wave compute times in this study . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.10 Coupled animation compute times in this study . . . . . . . . . . . . . . . . . . 59

7
 LIST OF FIGURES

1.1 A screenshot of Assassins Creed Origins taken by Esteb477 . . . . . . . . . . . . 13

2.1 A screenshot from Sea of Thieves . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1 Plots of SPH Kernels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2 An example of a Uniform Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 An example of a home cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.4 A visualization of the artificial downward force . . . . . . . . . . . . . . . . . . . 45
3.5 Coupled Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Particles applied to a section of the surface . . . . . . . . . . . . . . . . . . . . . 52
3.7 Coupled splash animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.8 Breaking wave animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.9 Boat Wake Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.10 Pure SPH Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.11 Pure Wave Surface Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.12 Pure SPH splash animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.13 Coupled Splash Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

8
 LIST OF SYMBOLS

CP U Central Processing Unit

GP U Graphical Processing Unit
SP H Smoothed-particle Hydrodynamics
SSBO Shader Storage Buffer Object

9
 ABSTRACT

This thesis project addresses the need for an interactive, real-time water animation tech-
nique that can showcase visually convincing effects such as splashes and breaking waves while
being computationally inexpensive. Our method couples SPH and wave equations in a one-
way manner to simulate the behavior of water in real-time, leveraging OpenGL’s Compute
Shaders for interactive performance and a novel Uniform Grid implementation. Through a
review of related literature on real-time simulation methods of fluids, and water animation,
this thesis presents a feasible algorithm, animations to showcase interesting water effects,
and a comparison of computational costs between SPH, wave equations, and the coupled
approach. The program renders a water body with a planar surface and discrete particles.
This project aims to provide a solution that can meet the needs of various water animation
use-cases, such as games, and movies, by offering a computationally efficient technique that
can animate water to behave plausibly and showcase essential effects in real-time.

10
 1. INTRODUCTION

Water animation has been a challenging task and an active area of research for many
years, due to its popularity and frequent use in special effects for feature films, TV, and video
games (Seymour, 2011). A wide range of methods have been developed to simulate water,
including particle-based techniques, wave equation techniques, and hybrid approaches that
combine multiple methods. However, each method has its own limitations and strengths,
and no single method can capture all the complex phenomena that occur in water.
A simulation tries to be exact with water's physical properties. However, an animation
of water tries to approximate the physical properties in order to provide a visually convinc-
ing version of it. Dynamic simulation of fluids has visually convincing techniques as well
as physically accurate techniques. Each technique has pros and cons and work better for
different use cases, but they can be combined to take advantage of their individual pros while
also reforming their cons.
Smoothed Particle Hydrodynamics (SPH) and wave equations are two widely used meth-
ods for water simulation and animation. SPH has the advantage of being highly versatile and
capable of modeling a wide range of water behaviors, including splashes and fluid dynamics.
Wave equations, on the other hand, are based on mathematical models that describe the
physics of waves and can produce highly realistic water surfaces in real-time.
Uniform Grids are also an important algorithm to take into account since it is an efficient
method for real-time collision detection. The Uniform Grid algorithm assumes that all
objects are approximately the same size and the concept involves creating a consistent grid,
where each individual cell of the grid is designed to be equal to or larger than the largest
object. We assign each object to one cell of the grid (the cell that the object currently
occupies), and then look at the neighborhood for any overlap with nearby objects. When
running this algorithm in parallel, every thread executes roughly the same amount of work,
so the execution divergence is low. Each thread in a given vicinity is designed to access
similar or the same memory locations as its neighboring threads, which results in minimal
differences in data access patterns, commonly referred to as data divergence (Karras, 2012).
For a particle-based fluid simulation, this algorithm works well but, it may suffer in other

11
situation like the “teapot-in-a-stadium” problem. The problem arises when a small object,
such as a teapot, is placed inside a larger, more complex environment, such as a stadium.
Traditional collision detection algorithms, such as those based on bounding volumes, can fail
to detect collisions between the teapot and the stadium walls. This is because the teapot
can be entirely contained within the stadium, and its bounding volume can be smaller than
the stadium’s bounding volume, making collision detection difficult.
This project has developed a new approach that couples SPH, and wave equations using
a uniform grid to achieve interactive real-time water animations with the ability to produce
effects such as breaking waves, boat wakes, and splashes while being computationally inex-
pensive. In this approach, wave equations are used to compute the surface waves of the water
body, while SPH is used to model the volumetric fluid dynamics of the water body. This
coupling allows for the creation of visually convincing water animations that are interactive,
computationally efficient, and behaviorally plausible.

1.1 Problem

The challenge of animating behaviorally plausible water has been a longstanding problem
in computer graphics. While various techniques exist to simulate water and produce visually
convincing results, many of these methods are either computationally expensive or not flex-
ible enough to accommodate different use cases. The purpose of this thesis is to develop a
method that couples SPH and wave equations to create a plausible, visually convincing water
animation that can showcase important effects such as splashes, boat wakes, and breaking
waves in real-time.
Existing SPH-based methods, which have been widely used in computer graphics for fluid
simulation, are highly versatile but can be computationally expensive. On the other hand,
wave equation-based methods that use a surface mesh provide a smooth representation of the
physics of waves, but cannot model certain water behaviors, such as splashes, and breaking
waves due to mesh topology. By coupling these two methods, we can harness the advantages
of each while minimizing their limitations.

12
 The main challenge in achieving this goal is to ensure that the resulting animation is
both physically plausible and computationally efficient. Real-time performance is important
for interactive applications, such as video games and simulators.

1.2 Purpose

Figure 1.1. A screenshot of Assassins Creed Origins taken by Esteb477

In figure 1.1, we can see different elements of a typical game scene – terrain, water, flora,
fauna, AI, global illumination, clouds, and fog to name a few. Since so many elements make
up a single scene, it is important to make sure that no single element takes up a majority
of the available computational resources. And since water can make up a large portion
of a scene, it is important to make sure that the water animation is not computationally
expensive because this will negatively impact the other elements.
The purpose of this thesis is to propose a new approach to real-time water animation
that combines Smoothed Particle Hydrodynamics (SPH) and wave equations. This method

13
aims to provide a plausible and efficient way to animate important water effects, such as
breaking waves and splashes, while keeping the computational cost within reasonable limits.
Our approach is not intended to replace other methods but rather to provide an alterna-
tive solution for specific use-cases where interactive performance is a priority. The purpose
of this thesis is to evaluate the feasibility of the proposed method, demonstrate its capa-
bilities through sample animations, and compare its performance to other water animation
methods.

1.3 Significance

Although physically-based simulations of water using SPH can provide realistic results,
it usually takes a lot of time and needs a lot of memory due to the number of particles
needed (Price, 2011). The simulation becomes better with more particles but every time the
particle count is increased, the memory needed by the simulation goes up significantly. This
may not matter to someone running the simulation on an extremely powerful PC that has
the space, and computational power it needs, but it does matter to a portion of the general
audience who do not have the resources needed to use SPH with many particles.
Wave equations have been used to model water in real-time applications because of their
speed and simplicity. They are based on the physics of water waves and can accurately
simulate the motion of water surfaces in real-time. Wave equations also have the ability
to model reflections and refractions of waves, making them a popular choice in game de-
velopment and virtual reality simulations (Yu et al., 2019). The main advantages of wave
equations is their efficiency. Wave equations are computationally inexpensive and can be
used to simulate large water bodies in real-time. However, there are also several challenges
associated with the use of wave equations for water animation. One challenge is the difficulty
of modeling volumetric effects such as foam and splashes, which are important for creating
visually convincing water animations.
Uniform grid is a widely used method for collision detection and resolution in fluid sim-
ulations, particularly in water animation. The method is based on dividing the simulation
domain into a uniform grid of cells, where each cell is a fixed size that is larger than the

14
size of the fluid particles. By grouping the particles into the cells they occupy, the method
reduces the number of pairwise collision checks between particles, resulting in significant
performance improvements. This approach is particularly useful in large-scale simulations,
where the number of particles can be very high, and efficient collision detection becomes a
significant bottleneck. In this way, uniform grid offers an efficient and effective way to han-
dle particle collisions in fluid simulations, making it an essential tool for the development of
realistic water animations. This project features a novel Uniform grid implementation.

1.4 Research Questions

The following questions needs to be answered to develop the deliverable:

• How to couple SPH and wave equations for animation of water?

• What are the performance, and plausibility tradeoffs of this coupled technique as
compared to pure SPH and wave equation implementations?

• How to implement uniform grid using compute shaders?

1.5 Hypothesis

Water animation using coupled SPH and wave equations will be able to produce a real-
time, visually convincing, and behaviorally plausible animation of water.

1.6 Project Deliverables

The deliverables of this thesis include:

• Algorithm and technique for water animation that couples SPH and wave equations.

• Analysis of the algorithm (Computation time vs Number of Particles).

• Representative images and videos of the technique at work.

• Compute shader uniform grid implementation.

15
1.7 Limitations

The limitations of the project include:

• Computational capabilities of the system used for this research.

• Drawbacks of SPH and wave equations (eg. lack of vortices).

• SPH and wave equations are not full Navier-Stokes simulations of water.

1.8 Delimitations

The delimitations of the project include:

• Rendering is not an aim and will not be addressed. The water body is rendered as
spheres, and a tessellated plane is used for the wave equation.

• Only a one-way coupling of SPH and wave equations is considered in order to main-
tain a fast, high-resolution animation.

1.9 Definitions

Computational Fluid Dynamics (CFD) is a numerical simulation technique that

uses mathematical models to analyze and solve problems related to fluid flow and heat
transfer.

Smoothed-Particle Hydrodynamics (SPH) is a numerical simulation method for

fluid mechanics that uses a Lagrangian approach to discretize and solve equations of
motion based on the properties of individual particles in the fluid.

Navier-Stokes Equation are a set of partial differential equations that describe the
motion of fluid in terms of its velocity, pressure, and viscosity, and are fundamental to
the study of fluid dynamics.

API is a set of functions provided to programmers to aid development.

16
 OpenGL (Open Graphics Library) is an open-source, cross-platform graphics API
that enables developers to create high-performance 2D and 3D graphics applications
for desktop and mobile devices.

Wave Equations are mathematical models that describe the behavior of waves in
terms of their amplitude, frequency, wavelength, and propagation speed.

Data Structures the organization of data (and its storage allocations in a computer).

Simulation is the process of creating a computer-based model or replica of a real-world

system or phenomenon to study its behavior or performance under different conditions
or scenarios.

Heightfield (heightmap) are two-dimensional representations of a three-dimensional

terrain, where each point on the terrain is represented by a single value indicating its
height or elevation above a reference level.

1.10 Summary

Fluid simulation and water animation have been important areas of research for decades,
with significant progress made using both SPH and wave equations. However, each of these
approaches have limitations in accurately and efficiently animating water for different ap-
plications. The main problem addressed in this research is the need for a water animation
technique that balances computational feasibility with visual accuracy. A pure physically-
based simulation of water is often computationally expensive for many applications, while
less accurate techniques may not produce a visually convincing animation. Therefore, the
purpose of this study is to develop a new method that couples the strengths of SPH and
Wave Equations to produce an efficient, visually pleasing, and behaviorally accurate anima-
tion of water. By combining these two methods, we aim to create a technique that is suitable
for a wide range of applications and can accurately simulate effects such as splashing, and
breaking waves in real-time.

17
 2. LITERATURE REVIEW

Water animation has been an intimidating task since the early days of computer graphics.
It is widely used in games, movies, and scientific simulators; sometimes making up a large
portion of the environment. The problem of water animation has been enduringly relevant
with numerous advances that have been made but there is still a long way to go. Initially,
water bodies were used as decorative pieces but they are now interactive areas that are
crucial to the application and its environment. Games like Assassin’s Creed: Origins, and
Sea of Thieves were recognized for the quality of their water and how it responds to the user
interacting with it (Gilroy, 2018). However, this quality is not achievable by all users as it
is heavily dependant on the power of their computer. Visual accuracy, and interactivity are
two things water animation often tries to balance; but one is usually sacrificed for the other.
The main philosophy driving most research in Water Animation is the lack of a “single model
fitting all needs” (Jensen & Golias, 2001).

Figure 2.1. A screenshot from Sea of Thieves

18
 The challenge lies in striking a balance between the computational complexity of phys-
ically realistic simulations, and the visual fidelity of animations that respond effectively to
user interaction. The need for accurate and visually convincing water animation is critical in
a wide range of fields, including games, movies, biological sciences, physics, and astronomy
(Karpinskyj & Sousa, 2013; Müller-Fischer et al., 2011; Zou et al., 2016).
In the past, researchers have explored various techniques for simulating and animating
water, including SPH and wave equations. While these methods offer distinct advantages,
they are not without limitations. For instance, SPH can produce physically accurate simu-
lations but requires large amounts of memory and computational resources as the number of
particles needed to simulate water increases (Price, 2011). On the other hand, wave equa-
tions provide a fast method for animating water but cannot accurately model volumetric
effects such as foam, splashes, and breaking waves (Yu et al., 2019).
To address these limitations, recent research has focused on combining SPH and wave
equations to produce fast, visually convincing, and physically accurate water animations. For
instance, Chentanez et al., (2015) presents a novel technique for simulating large scale liquid
phenomena in real time. The proposed method uses a combination of Eulerian, heightfield,
and particle methods to achieve accurate and efficient simulations of large bodies of water.
However, further work is needed to develop a more efficient and effective method for coupling
these techniques. This thesis study aims to contribute to this ongoing research by extending
current knowledge in the field and building on top of existing techniques to produce a more
effective method for water animation.

2.1 Water Animation using SPH

One of the main advantages of SPH is its ability to handle complex fluid behaviors, such
as splashes, foam, and turbulence. This is because SPH can simulate the interactions between
individual particles, which allows it to capture the details of fluid motion accurately. SPH
is also well-suited for parallel computing, which means it can be used to simulate large-scale
fluid systems efficiently. Becker, & Teschner (2007) mention that simulating breaking waves

19
and splashes are among the most desired features in fluid animation and that Lagrangian
methods like SPH are a promising way to capture these properties.
However, there are also several drawbacks and challenges associated with the SPH method.
One issue is the computational cost of SPH simulations can be high due to the large number
of particles that need to be simulated. Another challenge is the difficulty of handling bound-
ary conditions, such as fluid-solid interactions, which can be computationally expensive and
challenging to implement correctly.
To address some of these challenges, researchers have proposed various improvements
and extensions to the SPH method. For example, some studies have explored the use of
particle-mesh methods, which can reduce the computational cost of SPH simulations by
approximating the fluid field using a grid. Other studies have proposed hybrid methods
that combine SPH with other techniques, such as vortex methods or grid-based methods, to
improve the accuracy and efficiency of fluid simulations.
Lee & Han (2010) introduce “a 2D particle-based approach to achieve realistic water
surface behaviors for interactive applications” (p. 865). They argue that most physically
accurate simulations are slow because they need to compute Navier-Stokes equations, and
computational fluid dynamics (CFD) which are both time consuming for a large expanse
of water. Simplifying complex calculations to maintain interactive frame rates is a common
thread across most simulations of this kind. The authors split the water body into two parts:
the surface is simulated using wave equations and three-dimensional techniques whereas the
volume is simulated using a two-dimensional technique. Using a two-dimensional technique
for one area, and a three-dimensional technique for another consumes less resources than
using a three-dimensional technique for the entire simulation. However, there is no surface
for their water to flow on; their water body is “hanging in air”, this reduces the accuracy of
the interactions since free-floating water bodies do not occur in nature.
Along the same lines as Lee & Han (2010), Wang et al., (2021) mention their motivation
as the hindrance of large numerical integrations. Wang et al., (2017) do not try to optimize
an existing physically-based simulation method or make a new one, they instead use video
reconstruction. It is interesting to note that although they have the same motivation as
Lee and Han (2010), their approaches to solving the problem are vastly different from each

20
other. From Wang et al., (2017), Nie et al., (2021), and Nie et al., (2020) we can say that
video reconstruction is not widely used but is also not an entirely uncommon method for
fluid simulation; combining it with SPH has also been explored by Nie et al., (2020). This
method gives satisfying results, but it is not as widely usable as the more popular physically-
based methods as they require a source video for the simulation. Additionally, this model
entirely depends on the original reconstructed height field (Wang et al., 2010) and even a
small mistake or artifact could cause disastrous results. Pairing that with the fact that this
method cannot be rendered in real-time makes it impossible to use in interactive applications.
Moving away from animation methods and into attributes of water, we have capillary
waves ((a type of wave that occurs on the surface of a liquid where the restoring force is the
surface tension between the liquid and the surrounding air or another liquid), which are one
of the most common phenomena in fluid interaction but are notoriously “difficult to simulate
due to their fast travelling speed and high frequency” (Yang et al., 2016, p. 29). Enriching a
meshless technique like SPH with wave details has also been explored by Jang et al., (2010),
and Yang et al., (2016) but incorporating dynamic waves accurately into a low-frequency
result is difficult as the wave motion is simulated on its own (Yang et al., 2016), another issue
that pops up is that while gravity waves (waves caused by the oscillation of fluid particles
due to the restoring force of gravity and propagate through the entire depth of the fluid)
are handled effectively, capillary waves are not, since they arise in different situations and
have different wave equations (Yang et al., 2016). This brings to light the intricacies and
complexities of different waves in fluids that need to be simulated accurately in order to have
a plausible wave animation. Yang et al., (2016) mentions that they are simulating capillary
waves by using surface compression waves. Their assumption that it is a good approximation
is only considered reasonable because of the results of their experiment. The drawback here
is similar to Lee and Han (2010) where there are no boundary conditions, this limits the
accuracy and effectiveness of the model.

21
2.2 Water Animation using Wave Equations

Water animation has been an active research area in computer graphics, and the use
of wave equations has been a popular approach for simulating waves. The wave equation
method involves computing the velocity and height of waves at each point in a domain, and
then using this information to generate the motion of the water surface.
Jeschke and Wojtan (2015) present a new method for animating water waves by inter-
polating parameters across wavefronts, which allows for the creation of a wide variety of
wave shapes and behaviors. The method is based on the observation that the properties of
wavefronts, such as their amplitude and direction, change gradually across the surface of the
water. By using this property, the authors propose a wavefront parameterization technique,
which can be used to generate an animation of the entire water surface by interpolating
parameters between neighboring wavefronts. This approach enables the generation of com-
plex, and realistic water simulations with a high degree of control over the wave shapes and
behaviors. However, Jeschke and Wojtan note that this approach may not be suitable for all
types of water simulations, particularly those that require a high degree of physical accuracy.
Despite these challenges, the wave equation method has several advantages. It can pro-
duce visually appealing results and is well-suited for simulating large bodies of water, such
as oceans. Additionally, the spectral method allows for efficient simulation of waves with
different frequencies and amplitudes (Ducrozet et al., 2016).
Simulating water waves for a large expanse of water is computationally expensive but
the need for it grows each day as movies and games make larger environments that are com-
pletely filled with water. Schreck et al. (2019) “aims to simulate detailed water ripples down
to the level of individual raindrops and pebbles on a shoreline in an infinite, open ocean“
(p. 130:1); and introduce a novel approach for water wave animation using fundamental
solutions. The fundamental solutions, which are analytical solutions of the wave equation,
are precomputed and used to compute the water surface height and velocity in real-time.
The authors presented a technique to construct a library of fundamental solutions that is
both efficient to compute, and of high quality. The method is based on representing each
fundamental solution using a linear combination of spherical harmonics, and the coefficients

22
of the linear combination are computed using a fast Fourier transform. The authors demon-
strated the effectiveness of their method by simulating a variety of water phenomena such
as waves, wakes, and splashes. The main contribution of this work is the introduction of
a new approach for simulating water waves that is computationally efficient and produces
visually pleasing results. The method has several advantages over traditional approaches,
including the ability to handle complex geometries and the ability to simulate waves at
high resolutions. However, the method has some disadvantages, such as the requirement for
precomputation of the fundamental solutions, which can be time-consuming for large-scale
simulations. In addition, the method is currently limited to small-amplitude waves.

2.3 Other Methods of Water Animation

In addition to SPH and wave equations, other methods have been proposed for simulating
and animating water in computer graphics. These methods often aim to strike a balance
between physical accuracy, visual realism, and computational efficiency. In this section, we
review some of the other methods that have been proposed for water animation. These
include mesh-based techniques, and particle-based methods other than SPH. We present
their contributions, and discuss their respective strengths and weaknesses.

2.3.1 Ships, Splashes, and Waves on a Vast Ocean

Similar to pure SPH and wave equation approaches to water animation, coupled ap-
proaches have also been an area of active research for quite some time. These approaches
couple a particle method similar to SPH (or SPH itself) and a mesh method similar to wave
equations in order to produce realistic simulations of water.
Huang et al., (2021) proposes a new method for real-time simulation of water with highly
detailed breaking waves, splashes, and interactions with floating objects such as ships. They
couple Fluid-Implicit Particle (FLIP), and Boundary Element Method (BEM) to produce
realistic oceans by hybridizing “a local high-resolution FLIP solver and a surface-based BEM
solver for multiscale ocean modeling” (p.203:12). As the title describes, this method is
suitable for open water simulation (like vast oceans) with both water splashes near the main

23
object (like a boat), and waves that travel long distances. The authors test their method on
a large-scale ocean scene with multiple ships and waves. They compare their method with
previous approaches, such as using SPH alone or using a combination of SPH and grid-based
methods, and show that their method produces more realistic and detailed results while
maintaining real-time performance.

2.3.2 Water Surface Wavelets

Jeschke et al.’s (2018) paper “Water Surface Wavelets” proposes a novel method for
simulating water surfaces using wavelets. The authors’ proposed method uses a hierarchical
wavelet basis to represent the water surface, which allows for more efficient computation
and the ability to simulate small-scale features of waves. They also introduce a new wavelet
basis that is specifically designed to represent water surfaces, allowing for more accurate
simulation results.
The paper presents several examples of their method in action, including simulating water
droplets on a leaf, a ship sailing on a calm sea, and a waterfall. The results show that their
method is capable of producing realistic and visually appealing water simulations, while also
maintaining computational efficiency.

2.3.3 Real-Time Cartoon Water Animation

Most of the animation techniques here, and most fluid simulation techniques being de-
veloped focus on realism. Yu et al. (2007) however, deals with cartoonish water animation.
Computer Graphics research currently utilizes computational fluid dynamics for physically-
based simulations and while it produces good results, it cannot include an art of abstraction
such as cartoon stylization easily (Yu et al., 2007, p. 405). It is for this art of abstraction,
and the high cost of having professional animators accomplish this task that the authors
work to give the water stylized dynamics and representation. The drawback of this paper is
that it is template-based i.e., it can only be applied to classified types and not to a generic
water body.

24
2.4 Summary

Water animation is a topic that is constantly being worked on with numerous publications
and developments around it on a regular basis. With a wide variety of methods that can be
used and no method being singularly superior, often what sets a model apart is its use-case
and application. A model that needs to be applied for a shallow fountain has different needs
than one that needs to be applied to an ocean (although their behavior might be similar).
While Lee and Han (2010), Wang et al. (2017), and Yang et al. (2016) incorporate SPH to
animate water, Schreck et al. (2019) use wave equations, and Huang et al. (2021), Jeschke
et al. (2018), and Yu et al. (2007) use completely different approaches that often combines
multiple proven techniques.
A common thread across all these works are the authors’ need to solve the Navier-Stokes
Equation for water simulation if they use SPH. Most of the work here, and most of the
work done in the realm of modelling water phenomenon reduce the fluid into a surface repre-
sentation, however “they are unable to easily deal with complex three dimensional behaviors
such as flow around objects and dynamically changing boundaries” (Enright et al., 2002,
p. 737). Many researchers turned towards using a two-dimensional approximation of the
three-dimensional Navier-Stokes equation to make a dynamic fluid environment. While this
approximation is cost- and time-effective and also gives us a plausible dynamic simulation, a
lot of important data that could be used to make the simulation better is either disregarded
or lost.
By surveying the different approaches and their respective strengths and weaknesses, we
hope to have provided a comprehensive overview of the state of the art in water animation.

25
 3. METHODOLOGY

This methodology chapter presents the framework of the thesis project. It includes
workflows, and implementation details of our water animation model. SPH, and Wave
equations are discussed in some detail and our approach to coupling them is also explained.
Possibilities for each workflow, implementations, and animations are explored and discussed.

3.1 Development Tools

This project is implemented using C++, and OpenGL version 4.6. It also makes use of
OpenGL’s Compute Shaders. Installation of the corresponding OpenGL version and having
an Nvidia GPU card are required to run this project.

3.2 Research Approach

This research is developmental. A water animation is created by developing, exploring dif-

ferent implementations, and findings based on the methodologies discussed in the Literature
Review. The algorithmic implementation makes use of Compute Shaders to take advantage
of parallel processing in order to implement the SPH and wave equation calculations.
To show the visual and computational differences of the coupled approach, this project
has also implemented a pure SPH animation, and a pure wave equation animation. The
computational cost of the pure animations are then compared with the coupled implemen-
tation. The coupling used here is one–way — the wave surface affects the SPH particles but
the particles do not affect the surface.
The animation currently has a body of water made up of a flat planar surface with many
particles on it. When the user starts the animation, a wave is started from a point on the
surface and this in turn affects the particles and causes effects like splashes and breaking
waves. Parameters like number of particles, wave speed, viscosity, wave aatenuation etc., are
modifiable by the user via ImGui widgets so that they can fine tune the animation to their
liking.

26
3.3 SPH

First introduced by Gingold and Monaghan (1977), SPH is a CFD method that does not
require a mesh structure, and instead utilizes individual particles to represent the fluid flow
within a given system. The method employs a Lagrangian approach to model the behavior
of these particles as they move within the fluid.. Although it was originally meant for
modelling astrophysics phenomena it has now been extended for modeling fluid phenomena.
The method partitions the fluid volume into discrete particles. The boundaries of a specific
particle are blurred in order to create an even distribution of physical properties that are
associated with it (Kim, 2017). The main equation behind SPH is the integral approximation
of a field function given by (Monaghan, 2005)
Z
A(r) = A(r0 )W (r − r0 , h)dr0 (3.1)
ω

where r is a point in space, ω is the volume that contains r, and W is a kernel function with
a smoothing radius h.
The integral representation can be transformed into a particle approximation by a process
called discretization. The discretization involves a summation of all particles and is described
by the following
X
A(~r) = Aj Vj W (~r − r~j , h) (3.2)
j

where j loops over all particles, Aj is the field value at coordinate rj , and Vj is the volume
of particle j. Volume is the mass of an object divided by its denisty

m
V = (3.3)
ρ

Using this in equation 3.2, we get

X m
A(~r) = Aj W (~r − r~j , h) (3.4)
j
ρ

27
3.3.1 Density

SPH approximations can be utilized to solve the Navier-Stokes equations in a Lagrangian

formulation, which can be achieved through a two-step process. The continuity equation
can be omitted if the mass is fixed for each particle, and the density of each particle is
approximated with equation 3.4. This gives us

X m
A(~r) = ρj W (~r − r~j , h) (3.5)
j
ρ

which can then be simplified to

X
A(~r) = mj W (~r − r~j , h) (3.6)
j

3.3.2 Pressure

Due to the computational complexity of solving the Poisson pressure equation for near
incompressibility, it can be a time-consuming process. We overcome this by approximating
the pressure Pj of particle i with the equation of state to give weak compressibility

P = RρT (3.7)

where P is pressure, R is the molar gas constant, ρ is density and T is temperature. The
equation can also be written as
P = kρ (3.8)

where k is a gas stiffness constant that is dependant on T . This equation can be modified
as per the suggestion of Desburn (1996)

Pi = k(ρi − ρ0 ) (3.9)

where ρi is the density of particle i and ρ0 is the resting density.

28
3.3.3 Acceleration

The acceleration of a fluid particle can be expressed as

D~v ∂~v ∂~v ∂~v ∂~v ∂~v

~a = = + ~v · ∇~v = +u + +w (3.10)
Dt ∂t ∂t ∂x ∂y ∂z

however, since the particles move with the fluid, the convection term ~v · ∇~v is not needed.
This gives us
d~v F~i
a~i = = (3.11)
dt ρi

where v~i is the velocity of the particle i, F~i is the force density field of particle i, and ρi is
the density field at the location of particle i.

3.3.4 Forces

The momentum equation is the Navier–Stokes equation for incompressible, isothermal

Newtonian fluids with constant velocity. This is given by

∂~v
ρ( + ~v · ∇~v ) = −∇p + ρ~g + µ∇2~v (3.12)
∂t

where ρ is the density, ~v is the flow velocity, p is the pressure, ~g is an external force density
field (including gravity), and µ is the viscosity of the fluid. Here the acceleration terms are
arranged on the left side of the equation and the force terms are on the right side. This can
be reduced to
D~v
ρ = −∇p + ρ~g + µ∇2~v (3.13)
Dt
and the continuity equation
~v · ∇~v = 0 (3.14)

Equation 3.12 can be written as

d~
vi pressure viscosity external
= F~i + F~i + F~i (3.15)
∂t

29
 pressure
The pressure force F~i can be calculated using equation 3.6 for the pressure term −∇p

pressure X Pj
F~i = −∇p(~
ri ) = − mj ∇W (~
ri − v~j , h) (3.16)
j
ρj

but this does not produce a symmetric force (action is 6= reaction). The kernel gradient
has a value of zero at its central point, and when calculating the pressure force between
particles, each particle only considers the pressure of its counterpart. Due to differences in
the pressures at each particle’s location, the pressure forces will not be symmetrical. To
overcome this issue, a symmetric approximation of the gradient, as described by Muller et
al. (2003), is used to calculate the pressure forces:

pressure X P i + Pj
F~i =− mj ∇W (~
ri − v~j , h) (3.17)
j
2ρj

This gives us the pressure force. Similarly we can calculate a symmetric viscosity force by
using velocity differences and equation 3.6

viscosity X vj + vi 2
F~i =µ mj ∇ W (~
ri − v~j , h) (3.18)
j
2ρj

3.3.5 Kernel Functions

(a) Poly6 (b) Spiky (c) Viscosity

Figure 3.1. Plots of SPH Kernels

The poly6 kernel is a sixth-degree polynomial kernel that is used to compute the density
of each particle. It has the advantage of being symmetric and positive definite, but it is

30
prone to clustering particles together under high pressure (Muller, et al., 2003). To solve
this issue, the spiky kernel was introduced (Desbrun & Gascuel, 1996). It is a first-degree
polynomial kernel that produces a gradient that is proportional to the distance between
particles, allowing for a better representation of the fluid’s surface.
The spiky kernel is used to compute the pressure forces between particles. When two
particles are in close proximity, their repulsive forces vanish, which can lead to particle
clustering. The spiky kernel solves this problem by producing a repulsive force that increases
with particle separation. The kernel has a compact support, meaning that its influence is
limited to a finite range.
Finally, the viscosity kernel is used to compute the viscosity forces between particles. It
produces a force that opposes the relative motion of particles, helping to smooth out the fluid
and prevent it from appearing too turbulent. The viscosity kernel is based on the Laplacian
operator, which measures the rate of change of the velocity gradient. It has a positive definite
Laplacian, ensuring that it produces forces that oppose relative motion between particles.

3.3.6 Neighborhood

Efficiently identifying neighboring particles within a given radius is crucial for SPH simu-
lations. This study uses the Uniform Grid approach and searches over all particles. This has
a computational complexity of O(n2 ). The time and space complexity of using a Uniform
Grid in particle collision simulations depend on the size of the grid, the number of particles,
and the resolution of the simulation. A breakdown of the time complexity is:

• The time complexity of building the Uniform Grid is O(n), where n is the number
of particles. This is because each particle must be assigned to a cell in the grid,
which takes O(1) time per particle.

• The time complexity of calculating interactions between neighboring particles in the

same cell is O(k 2 ), where k is the average number of particles per cell. This is
because every pair of particles in the same cell must be considered for interaction.

31
 • The time complexity of calculating interactions between particles in neighboring
cells is O(k 2 ∗ m), where m is the number of neighboring cells that must be checked.
This is because for each neighboring cell, every pair of particles from the current
cell and the neighboring cell must be considered for interaction.

Overall, the time complexity of using a Uniform Grid in particle collision simulations is
O(n + k 2 ∗ m).
The space complexity of a Uniform Grid is O(N ), where N is the number of cells in the
grid. This is because each cell in the grid requires a fixed amount of memory to store particle
data.
Therefore, using a Uniform Grid in particle collision simulations provides an efficient way
to calculate interactions between particles, but the time and space complexity depend on
the size of the grid and the number of particles. We use a grid size of 32 × 32 × 32 cells
which was found to be the most optimal for our project.

3.4 Wave Equations

The linearized shallow water equation (Layton & van de Panne, 2002) assumes that the
speed of water varies slowly in space:

∂u ∂h ∂v ∂h ∂h ∂u ∂v
+g = 0, +g = 0, + d( + )=0 (3.19)
∂t ∂x ∂t ∂x ∂t ∂x ∂y

Where u, and v are fluid velocities along the x, and y directions, h is the height of the water
surface, d is the depth of water, and g is the acceleration due to gravity. It is also supposed
that the water surface is parallel to the xy coordinate plane.
Nishidate and Nikishkov (2005) differentiate the first equation in equation 3.19 with
respect to x, second equation in equation 3.19 with respect to y, third equation in equation
3.19 with respect to t and then combine all three equations in order to obtain the following
second order differential equation containing only the height h:

∂ 2h ∂ 2h ∂ 2h
 
= gd + (3.20)
∂t2 ∂x2 ∂y 2

32
d being a constant in the above equation makes it coincide with with the two-dimensional
wave equation. However, we can tell from our everyday experience that this is not sufficient
for realistic water movement modelling since damping is not present. Lets take a boat moving
on water for instance. As the boat moves forward, the wakes that spread behind it decrease
with time. This can be modelled if we add damping to our wave equation (Nishidate, &
Nikishkov, 2005):
∂ 2h ∂ 2h ∂ 2h
 
∂h
2
− k = c2 + (3.21)
∂t ∂t ∂x2 ∂y 2
where c is the speed of the wave and k is the damping constant.
In line with the above equations, we use the two-dimensional wave equation, which
describes the behavior of waves propagating through a medium in two dimensions, given by:

∂ 2w ∂ 2w ∂ 2w
 
= c2 + (3.22)
∂t2 ∂x2 ∂y 2

where w represents the wave height at a given point in space and time, t represents time, x
and y represent the spatial dimensions, and c represents the wave speed.

3.5 Program Flow

Given the nature of the method, SPH has good scalability and is inherently parallelized;
these two characteristics make Graphical Processing Units (GPU), and Compute Shaders a
good way to implement them. In this project, each particle has the following properties— po-
sition, velocity, force, density, and pressure. Each property is used to calculate the functions
described above. Compute Shaders are used for ease of coupling.
A Compute Shader is a type of shader stage that is primarily used for general-purpose
computations, rather than rendering graphics. Although it is capable of rendering graphics,
its main purpose is to perform tasks that are unrelated to drawing triangles or pixels, such as
physics simulations, and data processing. This animation is developed using four Compute
Shaders — one to calculate the forces acting on the particles, one to calculate the density of
the particles, one to calculate particle’s positions and velocities from its force and check for
boundaries, and one for computing and evolving waves using the wave equation.

33
 The SPH solver in this project takes advantage of the height of the wave from the wave
equation and integrates it such that water effects like splashes, and breaking waves can be
animated with ease. This coupling of equations tries to approximate the properties of water
in a visually convincing manner.
The program starts with the preparation and initialization of the water body which is
made up of particles and a displaced, tessellated plane. The water body is still at first
with the particles undergoing a “dam break” above it. This scene preparation includes
allocating space for the particles, wave surface, textures, and calculation matrices. During
initialization, the resting values of the particles, and wave surface are computed. Different
parameters (like particle density, viscosity, smoothing radius etc.,) are modifiable before and
after the animation starts. Waves, splashes, and wakes can be started using buttons available
as ImGui widgets.

3.5.1 Data Management

The project implements a data structure for the particles. This data structure contains
all the necessary attributes for an SPH particle— position, velocity, acceleration, radius,
density, and pressure. The mass, radius, and resting density of the particles are defined as
constants since they are the same for all particles.
Resting density is the density of the fluid at rest, i.e., when there are no external forces
acting on it. The resting density is typically set as the target density for the fluid. During
the simulation, if the density of a particle deviates from the target density, a pressure force
is applied to the particle to restore its density to the target value. This ensures that the fluid
remains incompressible, meaning that its density remains constant even under deformation.
The resting density is an important parameter in SPH simulations as it affects the behavior
of the fluid. For instance, a higher resting density leads to a more viscous fluid, while a lower
resting density results in a less viscous fluid.
A given data structure is considered flat if its individual elements are stored in a sequen-
tial and contiguous block of memory. The use of such a data structure provides us with
the advantage of cache locality, which makes them efficient and fast to access and process,

34
especially when compared to other data structures that are not stored in a contiguous mem-
ory block. This feature of flat data structures allows for easier and quicker traversal of the
elements during data processing. (Gunadi, & Yugopuspito, 2018). Therefore, we use a flat
data structure for the particles like:

struct Particle
{
vec4 pos; // xyz = pos, w = radius
vec4 vel; // xyz = vel, w = pressure
vec4 acc; // xyz = force, w = density
};

A single particle of this structure will occupy 4 × 4 × 3 = 48 bytes since a vec4 consists of
four floating–point values (a single float is 4 bytes) and the structure comprises of three such
vec4 values. The following table shows the total memory needed for different particle counts
when using this structure:

Table 3.1. Memory needed for different particle counts

Particle Count Size (bytes)
4096 196,608
8192 393,216
16,384 786,432
65,536 3,145,728

Textures provide a way to efficiently store and access large amounts of data, such as
heightfield maps, they can also be easily loaded and manipulated using various OpenGL
functions, allowing for dynamic updates and real-time rendering. For the wave, we use a
texture as a heightfield with the red channel of a pixel describing the height of the wave at
that point on the surface. A heightfield is a 2D grid of height values that represents a 3D
surface, with the height values defining the shape of the surface at each point. In this project,
the heightfield is used to simulate waves and other surface deformations by animating the
height values over time.

35
Compute Shader

From the official documentation (Compute Shaders, 2019), we know that compared to
other shaders stages like vertex and fragment shaders, compute shaders operate differently.
While other shader stages have well-defined built-in, and user-defined input values, compute
shaders do not. The rate at which a shader stage executes is determined by the type of shader
and the graphics hardware on which it is running, as well as the specific instructions and data
being processed by the shader; for instance, vertex shaders execute once per vertex, fragment
shaders execute once per fragment, and compute shaders execute once per invocation.
The “space” that compute shaders operate on is also known as the “thread group”; it is
determined by the number of thread groups and the size of each thread group specified when
the shader is launched. The total number of threads running in the shader is the product of
the number of thread groups and the size of each thread group, and each thread is assigned
a unique identifier that can be used to perform specific computations or operations. The
specific layout and organization of thread groups can also be customized by the programmer
using specialized functions and directives. Compute shaders lack user-defined inputs and
outputs, and the built-in inputs determine the location of a compute shader invocation in
the execution “space”.
Since there are no user-defined inputs and outputs for compute shaders; to get the nec-
essary data, compute shaders need to fetch it themselves, using textures, images, or shader
storage buffer objects (SSBO), which serve as the interface. Compute shaders must then
also write to one of these interfaces to perform calculations. In this project, an SSBO is used
for particle data, while a texture is used for the waves.

Work Groups

Within a compute shader, the work is partitioned into groups of threads that can execute
a specific number of threads concurrently. Each group is executed individually while the
threads within the group are processed in parallel. Similar to Gunadi and Yugopuspito
(2018), we parallelize the SPH algorithm by assigning a particle to each thread, and a given
thread calculates the SPH sums over that particle’s neighborhood.

36
 In OpenGL, work groups are automatically allocated to streaming multiprocessors (SMs),
and it is crucial to use the right number of work groups to achieve optimal performance. The
occupancy ratio, which is the number of threads that can run in parallel to the maximum
number of threads supported by the GPU, can have an impact on performance. Higher oc-
cupancy does not necessarily guarantee better performance, but low occupancyis guarenteed
to cause suboptimal performance. Therefore, care needs to be taken to choose the number
of work groups to achieve the best occupancy for their specific compute shader (Gunadi, &
Yugopuspito, 2018).

Shader Storage Buffer Objects

Shader Storage Buffer Objects (SSBOs) are an essential feature of modern graphics ren-
dering systems that allow data to be efficiently transferred between a compute shader and
other parts of the graphics pipeline (OpenGL, 2022). The main advantage of SSBOs is that
they provide a way to store a large amount of data in the graphics memory and to be accessed
by the shader programs. Their flexibility is particularly useful in simulation, and compu-
tational graphics applications where the number of data elements can change dynamically.
This makes SSBOs suitable for our project, which involves a large number of particles with
multiple attributes, and where the attribute values for each particle are expected to change
frequently.
SSBOs can also be written into by the shader program. This means that they can be
used to store the results of a computation, which can then be read by other parts of the
graphics pipeline. This is particularly useful in compute shaders, which can perform parallel
computations on a large number of data elements.
To use an SSBO, the data must first be written into a buffer object in system memory.
This buffer object can then be bound to the SSBO target, and the data can be accessed by
the shader program. The shader program can read from or write to the SSBO using the
built-in input and output variables provided by the GLSL. These input and output variables
define the data layout and can be used to access individual data elements within the SSBO
(OpenGL, 2022).

37
 Overall, SSBOs offer various ways to enhance the performance of graphics rendering
systems. In our project, we use SSBOs to implement the Uniform Grid, and the SPH
particle system. There are 3 SSBOs for the Uniform Grid, one for each array—count, start,
contents. For the SPH simulation, we use one SSBO that contains all the particles (and
their attributes) in our simulation. The inputs that are built-in to compute shaders offer a
convenient means of obtaining access to the data in the SSBOs.

3.6 Implementation

This project implements a novel compute shader approach to uniform grids, SPH, wave
equations. These latter two techniques are coupled in a one–way manner in order to produce
a visually plausible water animation. This section goes into the implementation details of
each approach and their coupling.

3.6.1 Uniform Grid

Figure 3.2. An example of a Uniform Grid

In a Uniform Grid, the simulation domain is divided into a regular grid of equal-sized
cells. Each cell is then assigned a unique index, which allows particles to be quickly located
and accessed during the simulation. When two particles are close enough to interact with
each other, their positions are checked to determine which cells they belong to. Only parti-
cles in neighboring cells are then considered for interaction, which significantly reduces the
computational cost of the simulation.

38
 Using a Uniform Grid for particle collision simulations offers several advantages. Firstly,
it reduces the computational complexity of the simulation by limiting the number of inter-
actions that need to be considered. This makes it possible to simulate larger numbers of
particles and longer timescales with a given amount of computational resources. Secondly,
it allows for efficient parallelization of the simulation, which can be split across multiple
processors or GPUs. This means that simulations can be run faster and more efficiently on
high-performance computing clusters.
However, there are also some limitations to using a Uniform Grid. In particular, the
grid size and resolution must be chosen carefully to ensure that particles are not missed
during the simulation. Additionally, the grid may not be well-suited to simulations where
the density of particles varies significantly across the simulation domain.
CPU Uniform grid implementations often use a 2D array of lists. Here, a cell is found
and the objects currently in that cell are inserted into that cell’s list (algorithm 2). When
trying to find which objects are in a cell, we first find a range of cells covered by the query
object and then return the corresponding overlapping objects (algorithm 3).
The following algorithms are examples of a Uniform Grid implementation on the CPU:

Algorithm 1 Grid Creation

1: N ← sizeworld /2
2: U nif ormGrid ← 2D list of size N × N

Algorithm 2 Grid Insertion

Require: U nif ormGrid is a 2D list of points in the Uniform Grid
1: cell ← CellCoord(point)
2: U nif ormGrid[cell.x][cell.y].push(point)

39
Algorithm 3 Check Overlap
Require: U nif ormGrid is a 2D list of points in the Uniform Grid
Require: aabb is an axis-aligned bounding box
Ensure: results is a list of points within the AABB
1: results ← empty list
2: for point p in U nif ormGrid[i][j] do
3: if aabb.contains(p) then
4: results.push(p)
5: end if
6: end for
7: return results

This approach has certain drawbacks to consider:

• Memory usage: The Uniform Grid requires a large amount of memory to store the
2D array of cells, and each cell may contain a list of objects that occupy that cell.
If the game world is large or the number of objects in the game is high, the memory
requirement can become prohibitive.

• Dynamic resizing: If the game world is dynamic, with objects being added or re-
moved, the Uniform Grid may need to be resized dynamically. This can be an
expensive operation, requiring the entire grid to be rebuilt and all objects to be
reinserted.

• Cache inefficiency: Accessing the 2D array of cells and the lists of objects within
each cell may not be cache-efficient, which can lead to slower performance.

An important point to note when implementing a Uniform Grid in the GPU is that we do
not have access to dynamic memory allocation. Therefore, we need to pre-allocate memory
and arrange for the data to be packed into memory. The drawbacks of the aforementioned
CPU approach can be circumvented by a GPU implementation.
Our novel compute shader implementation of a Uniform Grid uses four flat arrays that
work in tandem. We have an array each for count, start, and contents. The count array

40
describes how many particles are in each cell, and the start array describes what the starting
index of a particle is in the contents list. We make use of Parallel Prefix Sum to compute
start. Parallel prefix sum is a parallel algorithm used in computer science to efficiently
compute the prefix sum of an array of numbers. The prefix sum of an array is a new array
in which each element is the sum of all elements in the original array up to and including
that element.

Algorithm 4 AABB Increment Count

1: cellmin ← CellCoord(aabbmin ) . Get min coordinate of current AABB
2: cellmax ← CellCoord(aabbmax ) . Get max coordinate of current AABB
3: for i ← cellmin .x to cellmax .x do
4: for j ← cellmin .y to cellmax .y do
5: index ← Index(i, j)
6: atomicAdd(countindex , 1) . Increment count
7: end for
8: end for

Algorithm 5 Insert AABB into Uniform Grid cells

1: cellmin ← CellCoord(aabbmin ) . Get min coordinate of current AABB
2: cellmax ← CellCoord(aabbmax ) . Get max coordinate of current AABB
3: for i ← cellmin .x to clmax .x do
4: for j ← cellmin .y to cellmax .y do
5: index ← Index(i, j)
6: InsertAABB(aabb, index)
7: end for
8: end for

Axis-aligned bounding boxes (AABB) are made using the center of the given particle and
the smoothing length of the SPH simulation. The AABB increments the count (algorithm
4) for all cells the AABB overlaps.

41
Construction

• Compute count of objects in each cell. This is a simple sum of the number of objects
in a cell. Table 3.2 gives the count of each cell in figure 3.2.

• Compute start index of each cell. This is the exclusive parallel prefix sum (first
element of the array is excluded from the computation) of count and it can be
computed using a scan pattern.

Table 3.2. Cell index, count, and start arrays for grid in figure 3.2
Cell Count Start
0 2 0
1 2 2
2 0 2
3 3 4

• Insert objects into the contents list (using algorithm 6).

Table 3.3. content list for grid in figure 3.2

Index 0 1 2 3 4 5 6
Content 0 1 4 3 2 5 6

Algorithm 6 Insert objects

1: particle ← particlesi . Get current particle
2: cell ← CellCoord(particle) . Get the 2D cell coordinate of the current particle
3: index ← Index(cell) . Gets the 1D index of the cell
4: c ← atomicAdd(countindex , 1) . Increment count for this cell
5: s ← startindex . Get start value with this index
6: contentc+s = i . Insert point index

Table 3.2, and table 3.3 have been color coded as per the uniform grid cells in figure 3.2.

42
Collision Query

• For each grid cell (i, j) overlapped by current particle.

• Get count particles starting at start index.

• If particle in in overlapped area, proceed with collision resolution.

For a given object, we generate a list of the potentially colliding objects.

We can avoid self-collisions, and double collisions by generating a unique home cell for
each pair of objects. We only add collisions from the home cell for each pair of objects. We
define a home cell as the cell containing the maximum of the minima of the 2 AABBs. Using
either the minimax or maximin as the home cell will provide a unique cell to consider for a
given pair of objects.
It is important to note that our animation does not need a home cell since it deals with
points; which can only occupy a single cell at a given time. A home cell will be required
when dealing with rectangles or circles which can occupy multiple cells at a given time.

Figure 3.3. An example of a home cell

In figure 3.3, the home cell of AABB 0 and AABB 1 is cell 11.

43
Assumptions

We make the following assumptions for our Uniform Grid implementation:

• The object we are inserting is a given particle’s AABB.

• Multiple particles can be in a cell but a given particle can only be in one cell.

SSBO

We use three SSBOs, one for each of our arrays—count, start, and contents; all of them
contain a single integer array. The length of the count, and start arrays are independent of
particle count, their length will be the total number of cells in the Uniform Grid, and the
length of the contents array will be the total number of particles.
Since a single integer occupies 4 bytes, the aforementioned arrays will occupy the following
amount of memory for different particle, and grid cell counts:

Table 3.4. Memory occu- Table 3.5. Memory occupied

pied by contents array for by count, and start for different
different particle counts grid counts
Particle Count Size (bytes) Grid count Size (bytes)
4096 16,384 8×8×8 2048
8192 32,768 16 × 16 × 16 16,384
16,384 65,536 32 × 32 × 32 131,072
65,536 131,072 64 × 64 × 64 1,048,576

3.6.2 SPH

This project makes use a compute shader for its SPH solver. The compute shader serves
three purposes— calculate the density and pressure on the particles, calculate the forces
acting on the particles, and integrate everything in order to set the positions and velocities
of the particles as well as check boundary conditions.

44
Foam

Particles with low pressure are treated as foam particles since they have likely been
ejected from the water surface. These particles are also rendered with a lighter color to more
closely simulate the white of foam.

Force Computation

In addition to the previously described forces that act on the particles, we also exert an
artificial downward force on particles that are too high above the wave surface so that the
splashes do not fly too high, ensuring that the particles are contained on the wave.

Figure 3.4. A visualization of the artificial downward force

In a particle-based simulation like SPH, particles represent fluid elements that interact
with each other through forces. When a wave breaks, particles experience a large upward
force that can cause them to be ejected from the surface, leading to unrealistic splash be-
havior.
By applying an artificial downward force on the particles, we can counteract this upward
force and keep the particles close to the surface. This can help create a more realistic wave
behavior by preventing the particles from settling into the wave troughs as well as prevent
particles from flying too high and forming unrealistic splashes.

45
 However, it’s also worth noting that applying too strong of a downward force can also
lead to unrealistic behavior, such as waves that do not break properly or waves that appear
to “stick” to the surface.

Density and Pressure Computation

The density and pressure computation of the particles have a direct influence from the
wave in addition to the calculations previously mentioned. The density contribution from
the wave follows algorithm 10. The force contribution from the wave produces a pressure
and viscosity force that is computed using algorithm 8, and algorithm 11 similarly to the
density contribution.

Boundary Conditions

Many SPH implementations make use of clamped boundaries but this can have certain
drawbacks like:

• Artificial behavior: Clamped boundaries can create artificial forces and density gra-
dients that can affect the behavior of the fluid in the simulation. This can lead to
inaccurate results, especially in regions close to the boundaries.

• Discontinuities in pressure and velocity: When particles collide with the boundaries,
they can experience sudden changes in velocity and pressure, which can create dis-
continuities in the fluid flow. These discontinuities can cause instabilities in the
simulation and affect its accuracy.

We use periodic boundaries to overcome these challenges since they allow particles that
leave the simulation domain to re-enter the domain from the opposite side, as if the simulation
domain was repeated periodically in space. This helps to avoid the problem of particles being
artificially stopped at the boundary and causing inaccurate behavior.
Periodic boundaries can also help to avoid the formation of artificial waves at the edges
of the simulation domain, which can occur with clamped boundaries. When using clamped
boundaries, waves can reflect off the boundary and interact with incoming waves, leading to

46
a complex interference pattern. With periodic boundaries, waves can simply wrap around
the domain, avoiding this interference and resulting in a more natural simulation. This
is more useful since it lets us keep a constant number of particles on the wave. In our
implementation, the particles do not respawn exactly at the boundary, thereby preventing
thrashing (which can be a common issue in periodic boundaries).

Algorithm 7 Periodic Boundary

1: eps ← 0.05 . Small offset
2: if particlepos ≥ (boundarylower + eps) then
3: particlepos ← particlepos + (2 × boundarylower − 2 × eps) . Offset the position
4: end if
5: if particlepos ≤ (boundaryupper + eps) then
6: particlepos ← particlepos − (2 × boundaryupper − 2 × eps) . Offset the position
7: end if

Additionally, periodic boundaries can help to reduce the computational cost of a simula-
tion, since fewer particles need to be simulated to cover a larger domain.

SSBO

We set the maximum number of SPH particles for our project to be 64 × 1024 = 65536.
The particle SSBO is initialised to a buffer of these many particles so that it is easier to
update the number of particles at runtime. The size of a single particle is 48 bytes; therefore
the particle SSBO occupies 65536 × 48 = 31, 45, 728 bytes. This buffer will contain empty
elements if there are fewer than 65536 particles.

3.6.3 Wave Surface

The wave surface is a displaced, tessellated plane in this project. Each vertex of the
surface has its height calculated using the wave equation. The wave equation is computed
and written into an image whose red channel describes the height of the wave at that point.

47
This red channel value is scaled and mapped to the surface in order to show the wave moving
across it.
A compute shader is used to calculate the wave equation using finite-difference time-
domain (FDTD) method to discretize and implement equation 3.22. The simulation is
performed on a 2D grid represented by an image. The shader operates on three images
which represent the wave at times t, t-1, and t-2. We also define a “neighborhood” struct to
represent the color values of the current pixel, and the pixels immediately above, below, to
the left, and to the right of it.
The two main modes of operation are:

• Initialisation: Initializes the wave at time t with a simple wave or a splash.

• Evolution: Evolves the wave at time t based on the values of the waves at times t-1,
and t-2. The shader computes the new wave values using the FDTD method and
stores the result in an image.

The computation of the evolving wave at time t is based on a neighborhood of pixels

surrounding each pixel in the image and three constant attributes for the wave—Lambda,
Attenuation, and Beta. The shader first retrieves the colors of the six surrounding pixels
and then applies the FDTD method using the attributes to evolve wave values based on the
retrieved pixel colors. The result is stored in the output image.

48
3.6.4 Coupling

(a) At rest (b) Animated

Figure 3.5. Coupled Animation

The aforementioned SPH and wave equations are coupled by making use of the wave
heightfield in the SPH computation. Since the red channel of the image describes the height
of the wave at that point on the surface, it is straightforward to extract data from it for the
SPH computation.
Following algorithm 8, we first get the height of the wave at the particle’s position using
the particle’s X– and Z–coordinate.

Algorithm 8 Calculate wave height at a given particle position

1: function Wave Height(pos) . Particle position is passed in to the function
2: posxz ← posxz + waveshif t . Offset the position by the wave shift
3: posxz ← posxz × wavescale . Scale the position by the wave scale
4: uv ← 0.5 × posxz + vec2(0.5) . Calculate the texture coordinates
5: h ← texture(heightf ield, uv)r . Read the wave height from the texture
6: return (h − 0.5)
7: end function

49
 We make use of the wave normal (calculate using algorithm 9) to compute a wave force
and implement a one–way coupling where the wave affects the particles by exerting forces
on it. Therefore, even though the wave surface and the particles form a single water body,
they are treated as independent entities with the wave exerting external forces on the SPH
particles.

Algorithm 9 Calculate Wave Normal at particle position

1: function Wave Normal(pos)
2: posxz ← posxz + waveshif t . Offset the position by the wave shift
3: posxz ← posxz × wavescale . Scale the position by the wave scale
4: uv ← 0.5 × posxz + vec2(0.5) . Calculate the texture coordinates
5: duv ← vec3(1.0/texture_size, 0.0) . Calculate the texture coordinate offset
. Calculate the x component of the normal
6: dx ← texture(heightf ield, uv + duvxz )r − texture(heightf ield, uv − duvxz )r
. Calculate the y component of the normal
7: dy ← texture(heightf ield, uv + duvyz )r − texture(heightf ield, uv − duvyz )r
. Combine and normalize the x and y components of the normal
8: n ← normalize(vec3(dx, 0.1, dy))
9: return n . Return the wave normal
10: end function

Once we know the height of the wave at the particle’s position, we allow the wave to
exert pressure and viscosity forces on the particles. This is done in the SPH solver’s density
and pressure computation (algorithm 10), as well as in the force computation. We make
use of a cubic spline kernel (algorithm 11) to provide a smooth interpolation of properties
such as density and pressure from neighboring particles. The function is designed to be
differentiable, which allows it to be used to calculate forces such as pressure gradients.

50
Algorithm 10 Compute density contribution from wave
1: function Rho Wave(pos)
2: rho ← 0.0
3: h ← wave_height(pos) . Get wave height at this position
4: d ← posy − h
5: t ← (1.0, 0.0, 0.0) . Constant +X-Axis vector
6: b ← (0.0, 0.0, 1.0) . Constant +Z-Axis vector
7: if d < H then . Check if wave height is less than particle neighborhood radius
8: nb ← 5
9: for i = 0 to nb − 1 do
10: dx ← (i − nb/2) × particle_diameter
11: for j = 0 to nb − 1 do
12: dy ← (j − nb/2) × particle_diameter
13: b_pos ← p + dx × t + dy × b
14: b_posy ← wave_height(b_pos)
15: rho ← rho + 1.0 × W _cubic(p − b_pos)
16: end for
17: end for
18: end if
19: return rho
20: end function

51
Algorithm 11 Cubic kernel weight function
Require: r: distance between two particles
Require: H: smoothing radius
Require: kdim: kernel dimensionality constant
1: function W Cubic(r)
2: q ← r/H
3: if q ≥ 1.0 then
4: return 0.0
5: else if q ≤ 0.5 then
6: return kdim × (6.0 × q 2 × (q − 1.0) + 1.0)
7: else
8: q1 ← 1.0 − q
9: return kdim × 2.0 × q13
10: end if
11: end function

In our project, the SPH particles can be applied to certain parts of the wave surface to
focus the effects around that area (accomplished by shifting or scaling the surface). This
enables us to enrich sections of the wave with higher detail.

Figure 3.6. Particles applied to a section of the surface

52
Splashes

(a) Particle body dissected to see the wave

surface (b) Splash peak

Figure 3.7. Coupled splash animation

The camera is in a perspective view for the splash animation in order to see the splashes,
and foam. The wave starts with a depression in the middle of the wave surface which evolves
to animate a “splash”.

Breaking Waves

Figure 3.8. Breaking wave animation

53
 The camera has been moved to a slight orthographic view for the breaking wave animation
in order to see the wave form better. The wave seed is on an edge of the surface, the wave
then evolves towards the opposite edge.

Boat Wake

(a) Boat wake (b) With dissected particle body

Figure 3.9. Boat Wake Animation

The camera moves to a perspective view for the boat wake animation in order to see the
waves, wakes, and splashes. The wave seed is a single point at the boat’s starting location,
the wave then evolves from a straight line starting that that point. The wake form can be
seen by the white colored SPH particles, these trace the path of the wake.

3.7 Rendering

As mentioned before, the rendering of this animation is not an aim of this research but
care has been taken to ensure that the animation is rendered in a visually plausible manner
so that the visual properties of water are exhibited to some extent.
There are four pairs of vertex, and fragment shaders—one for the particles, one for the
wave surface, one for the boat, and one for the skybox,

54
3.7.1 Particles

(a) At rest (b) Animated

Figure 3.10. Pure SPH Animation

The particles are rendered as spheres using “ray-sphere intersection”, and their color
is computed using their density and pressure attributes. Due to lack of particles near the
surface the calculated pressure is low. This is indicated with a whiter particle color and can
be used for foam or spray rendering.

3.7.2 Wave Surface

(a) At rest (b) Animated

Figure 3.11. Pure Wave Surface Animation

55
 The wave surface is rendered as a flat diffuse plane using per-fragment Phong lighting.

3.8 Results

This section provides the visual animations, and computational time of our coupled
approached when compared to pure SPH and wave equation approaches.

3.8.1 Visuals

The following table shows a visual comparison between the different types of animation:

Table 3.6. Visual comparison of animations

SPH Wave Equation Coupled

Splash

Breaking Wave

Boat Wake

3.8.2 Performance

The performance of this project was measured on an Razer Blade 14 laptop using an
NVIDIA GeForce RTX 3060, and an AMD Ryzen 9 6900HX Processor.
The following tables show the time taken by various computations in our implementation:

56
 Table 3.7. Uniform Grid compute times in this study
Particle Count Time (ms)
8192 0.14
16384 0.14
32768 0.15
65536 0.15

A crucial factor in the computation cost of the Uniform Grid is the grid size. Table 3.7
had a grid size of 32x32x32, which was also found to be the most optimal size. At larger
grid sizes, clustering caused problems for particle-particle collisions.
With no Uniform Grid (or a grid size of 1x1x1), the time taken by the Uniform Grid
computation (with 8192 particles) is close to 3 ms whereas it is only around 0.18 ms when
the grid size is set to an optimal value. This makes the overall computation time taken by
the coupled animation increase from around 2 ms on average to around 6 ms, almost tripling
the computation time.

Table 3.8. SPH compute times in this study

Particle Count Time (ms)
8192 0.60
16384 1.10
32768 3.90
65536 5.30

57
 (a) With 8192 particles (b) With 16384 particles

Figure 3.12. Pure SPH splash animation

After considering the visual quality, and computational cost with various particle counts;
a particle count of 16384 was found to be the most optimal for the SPH animation. At
this value, the computational cost is low, the framerate is high, and all water effects can
be showcased effectively. Decreasing the particle count prevents the animation from show-
casing all our water effects effectively whereas increasing the particle count increases the
computational time with a severe decrease in framerate, and makes the animation extremely
unstable.

Table 3.9. Wave compute times in this study

Wave Resolution Time (ms)
256x256 0.04
512x512 0.06
1024x1024 0.23

The computation time of the wave is very low and almost negligible when compared to
the computation of the Uniform Grid, and the SPH animation.

58
 Table 3.10. Coupled animation compute times in this study
Particle Count Time (ms)
8192 0.50
16384 1.10
32768 1.80
65536 3.20

After experimentation considering visual quality, and computational cost with various
particle counts and wave resolutions; a wave resolution of 256x256, and a particle count of
8192 was found to be the most optimal for the coupled animation. At these values, the
computational cost is low, the framerate is high, and all water effects are shown effectively.
Decreasing the particle count prevents the animation from showcasing all our water effects
effectively whereas increasing the particle count increases the computational time with no
major change in the visual quality of the water effects, and makes the animation unstable.

(a) With 8192 particles (b) With 16384 particles

Figure 3.13. Coupled Splash Animation

59
3.8.3 Conclusion

Looking at the visuals, and the compute time breakdowns, we can make the following
conclusions:

• Our coupled approach is able to showcase water effects like foam, splashes, and
breaking waves using fewer particles, and in lesser time.

• A pure SPH approach is able to showcase water effects like splashes, and breaking
waves but requires at least double the amount of particles as the coupled approach.
Therefore, using up more memory, and having a higher computational cost.

• A pure wave equation approach using a surface mesh is incredibly fast but cannot
show any volumetric water effects.

60
 4. SUMMARY
4.1 Conclusion

The problem addressed by this research is that a pure SPH–based water animation is
computationally expensive, and a pure wave equation–based water animation does not accu-
rately simulate water behavior. By coupling SPH, and wave equation, we are able to develop
a water animation method that accurately approximates water behavior beyond surface rep-
resentations, and is computationally inexpensive. The need for such an animation increases
as the demand for fast, and realistic water animation for games, movies, and VFX increases
(Seymour, 2011). This research is unconcerned with reproducing water’s exact physical prop-
erties or have a hyper-realistic rendering of water. Instead, this research is concerned with
approximating water’s properties, and behavior to produce a visually convincing animation.
Our approach has several potential applications in areas such as video games, movies,
and VFX. We believe that our research provides valuable insights for the development of
more advanced water simulation techniques that can be used in a wide range of fields.
Our method makes use of Compute Shaders to aid in parellelization, implement a uniform
grid, solve wave equations, and implement SPH. The particles and wave surface are coupled
in a one-way manner where the wave surface affects the SPH particles. The deliverables of the
research are animations that shows water effects like foam, splashes, and breaking waves,
a novel compute shader uniform grid implementation, an algorithm for water animation,
and an analysis of the coupling. The computational cost of the coupled implementation is
compared with that of pure SPH and wave equation implementations, and is found to be
relatively cheaper. Our coupled approach is able to showcase water effects that are either
costly or impossible for pure SPH or wave surface implementations to showcase.

4.2 Future Work

Although the proposed method of coupling SPH, and wave equations for water anima-
tion is promising, there are several avenues for future research that can improve, and extend
the method. One direction is to focus on realistic rendering of the water surface. While

61
the current approach produces visually convincing results, further research can investigate
techniques to achieve physically accurate rendering that accounts for properties such as sur-
face tension, light transport, and caustics. Another area for future work is the simulation
of collisions between water and other objects, such as ships or rocks. This requires coupling
the water simulation with rigid body dynamics, which can pose challenges in terms of com-
putational efficiency and stability. Along with this, many demos are initialized with a “dam
break”, but that is not useful for games. Initial particle positions must respect the fluid rest
density, and must be initialized as such. Wrap-around effects due to the periodic boundary
condition of the SPH solver can also be prevented.

62
 5. REFERENCES

[1] Price, D. (2012). Smoothed Particle Hydrodynamics: Things I Wish My Mother

Taught Me. In Advances in Computational Astrophysics: Methods, Tools, and Out-
come.(pp. 249).

[2] M.S. Shadloo, G. Oger, & D. Le Touzé (2016). Smoothed particle hydrodynamics
method for fluid flows, towards industrial applications: Motivations, current state,
and challenges. Computers & Fluids, 136, 11-34.

[3] Gingold, R., & Monaghan, J. (1977). Smoothed particle hydrodynamics: theory and
application to non-spherical stars. Monthly Notices of the Royal Astronomical Society,
181(3), 375-389.

[4] Crespo AC, Dominguez JM, Barreiro A, Gómez-Gesteira M, Rogers BD (2011) GPUs,
a New Tool of Acceleration in CFD: Efficiency and Reliability on Smoothed Particle
Hydrodynamics Methods. PLoS ONE 6(6): e20685.

[5] Lee, H., & Han, S. (2010). Solving the shallow water equations using 2d sph particles
for interactive applications. The Visual Computer, 26(6), 865-872

[6] Wang, C., Wang, C., Qin, H., & Zhang, T. Y. (2017). Video-based fluid reconstruction
and its coupling with SPH simulation. The Visual Computer, 33(9), 1211-1224

[7] Yang, S., He, X., Wang, H., Li, S., Wang, G., Wu, E., & Zhou, K. (2016, July). En-
riching SPH simulation by approximate capillary waves. In Symposium on Computer
Animation (pp. 29-36).

[8] Schreck, C., Hafner, C., & Wojtan, C. (2019). Fundamental solutions for water wave
animation. ACM Transactions on Graphics (TOG), 38(4), 1-14.

[9] Wang, H., Miller, G., & Turk, G. (2007, August). Solving general shallow wave equa-
tions on surfaces. In Proceedings of the 2007 ACM SIGGRAPH/Eurographics sym-
posium on Computer animation (pp. 229-238).

63
[10] Yu, J., Jiang, X., Chen, H., & Yao, C. (2007). Realtime cartoon water animation.
Computer Animation and Virtual Worlds, 18(45), 405-414.

[11] Jensen, L. S., & Golias, R. (2001, September). Deep-water animation and rendering.
In Game Developers Conference (Gamasutra).

[12] Nie, X., Hu, Y., & Shen, X. (2020). Physics-preserving fluid reconstruction from
monocular video coupling with SFS and SPH.The Visual Computer,36(6), 1247-1257.

[13] Nie, X., Hu, Y., Su, Z., & Shen, X. (2021, September). Fluid Reconstruction and Edit-
ing from a Monocular Video based on the SPH Model with External Force Guidance.
In Computer Graphics Forum(Vol. 40, No. 6, pp. 62-76).

[14] An, H., & Park, J. (2015). A Lagrangian Approach on Video based Fluid Animation.
InThe First Workshop on Art, Culture, Game, Graphics, Broadcasting and Digital
Contents.

[15] Angst, R., Thuerey, N., Botsch, M., & Gross, M. (2008, October). Robust and efficient
wave simulations on deforming meshes. InComputer Graphics Forum(Vol. 27, No. 7,
pp. 1895-1900). Oxford, UK: Blackwell Publishing Ltd.

[16] Jang, T., Kim, H., Bae, J., Seo, J., & Noh, J. (2010). Multilevel vorticity confinement
for water turbulence simulation.The Visual Computer,26(6), 873-881.

[17] Gilroy, J. (2018, August 10). Why the Water in Sea of Thieves Is So Mesmerizing.
Fandom. Retrieved October 23, 2022, from https://www.fandom.com/articles/sea-of-
thieves-water

[18] Gingold, R. A., & Monaghan, J. J. (1977). Smoothed particle hydrodynamics: theory
and application to non-spherical stars. Monthly notices of the royal astronomical
society, 181(3), 375-389.

[19] Kim, D. (2017). Fluid engine development. CRC Press.

[20] Monaghan, J. J. (2005). Smoothed particle hydrodynamics. Reports on progress in

physics, 68(8), 1703.

64
[21] Desbrun, M., & Gascuel, M. P. (1996). Smoothed particles: A new paradigm for
animating highly deformable bodies. In Computer Animation and Simulation96 (pp.
61-76). Springer, Vienna.

[22] Müller, M., Charypar, D., & Gross, M. H. (2003, July). Particle-based fluid simulation
for interactive applications. In Symposium on Computer animation (Vol. 2).

[23] Gunadi, S. I., & Yugopuspito, P. (2018, August). Real-time gpu-based sph fluid simu-
lation using vulkan and opengl compute shaders. In 2018 4th International Conference
on Science and Technology (ICST) (pp. 1-6). IEEE.

[24] Obeid, N. M. (2011). Compact binning for parallel processing of limited-range func-
tions.

[25] Peachey, D. R. (1986). Modeling waves and surf. ACM Siggraph Computer Graphics,
20(4), 65-74.

[26] Chentanez, N., & Müller, M. (2010, July). Real-time Simulation of Large Bodies of
Water with Small Scale Details. In Symposium on Computer Animation (pp. 197-206).

[27] Jang, T., & Noh, J. (2015). A geometric approach to animating thin surface features
in smoothed particle hydrodynamics water. Computer Animation and Virtual Worlds,
26(2), 161-172.

[28] Fujisawa, M., Nakada, T., & Mikawa, M. (2017). Particle-based shallow water simula-
tion with splashes and breaking waves. Journal of information processing, 25, 486-493.

[29] Chentanez, N., Müller, M., & Kim, T. Y. (2015). Coupling 3D eulerian, heightfield
and particle methods for interactive simulation of large scale liquid phenomena. IEEE
transactions on visualization and computer graphics, 21(10), 1116-1128.

[30] Bailey, M. (2016). OpenGL Compute Shaders. Oregon State University.

[31] Sellers, G., Wright Jr, R. S., & Haemel, N. (2013). OpenGL superBible: comprehensive
tutorial and reference. Addison-Wesley.

65
[32] Compute Shader. (2019, April 22). OpenGL Wiki, . Retrieved 18:52, February 25,
2023 from http://www.khronos.org/opengl/wiki_opengl/index.php?title=Compute_Shader&ol-
did=14536.

[33] Ihmsen, M., Cornelis, J., Solenthaler, B., Horvath, C., & Teschner, M. (2013). Implicit
incompressible SPH. IEEE transactions on visualization and computer graphics, 20(3),
426-435.

[34] OpenGL. (2022). Shader storage buffer object. Retrieved February 26, 2023, from
https://www.khronos.org/opengl/wiki/Shader_Storage_Buffer_Object

[35] OpenGL. (2019). OpenGL 4.6 specification. Retrieved February 26, 2023, from
https://www.khronos.org/registry/OpenGL/specs/gl/glspec46.core.pdf

[36] Jones, D. (2014). OpenGL insights. CRC Press.

[37] Akenine-Mo, T., Haines, E., & Hoffman, N. (2018). Real-time rendering.

[38] Woo, M., Neider, J., Davis, T., & Shreiner, D. (1999). OpenGL programming guide:
the official guide to learning OpenGL, version 1.2. Addison-Wesley Longman Publish-
ing Co., Inc..

[39] Rost, R. J., Licea-Kane, B., Ginsburg, D., Kessenich, J., Lichtenbelt, B., Malan, H.,
& Weiblen, M. (2009). OpenGL shading language. Pearson Education.

[40] Tessendorf, J. (2001). Simulating ocean water. Simulating nature: realistic and inter-
active techniques. SIGGRAPH, 1(2), 5.

[41] Chan, K. H., Ke, W., & Im, S. K. (2018). Particlemesh coupling in the interaction of
fluid and deformable bodies with screen space refraction rendering. Computer Anima-
tion and Virtual Worlds, 29(1), e1787.

[42] Hinsinger, D., Neyret, F., & Cani, M. P. (2002, July). Interactive animation of ocean
waves. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on
Computer animation (pp. 161-166).

66
[43] Huang, L., Qu, Z., Tan, X., Zhang, X., Michels, D. L., & Jiang, C. (2021). Ships,
splashes, and waves on a vast ocean. ACM Transactions on Graphics (TOG), 40(6),
1-15.

[44] Jeschke, S., Skivan, T., Müller-Fischer, M., Chentanez, N., Macklin, M., & Wojtan, C.
(2018). Water surface wavelets. ACM Transactions on Graphics (TOG), 37(4), 1-13.

[45] Ducrozet, G., Bonnefoy, F., Le Touzé, D., & Ferrant, P. (2016). HOS-ocean: Open-
source solver for nonlinear waves in open ocean based on High-Order Spectral method.
Computer Physics Communications, 203, 245-254.

[46] Karras, T. (2012). Thinking parallel, part i: Collision detection on the gpu.

[47] Jeschke, S., & Wojtan, C. (2015). Water wave animation via wavefront parameter
interpolation. ACM Transactions on Graphics (TOG), 34(3), 1-14.

[48] Bridson, R. (2015). Fluid simulation for computer graphics. CRC press.

67