# Ships, Splashes, and Waves on a Vast Ocean

Libo Huang  KAUSTVisual Computing Center, Thuwal 23955, KSA ,  Ziyin Qu  UCLA
and UPenn ,  Xun Tan  Zenus TechnologyShenzhen ,  Xinxin Zhang  Zenus
TechnologyShenzhen ,  Dominik L. Michels  KAUSTVisual Computing Center, Thuwal
23955, KSA and  Chenfanfu Jiang  UCLA and UPenn

###### Abstract.

The simulation of large open water surface is challenging using a uniform
volumetric discretization of the Navier-Stokes equations. Simulating water
splashes near moving objects, which height field methods for water waves
cannot capture, necessitates high resolutions. Such simulations can be carried
out using the Fluid-Implicit-Particle (FLIP) method. However, the FLIP method
is not efficient for the long-lasting water waves that propagate to long
distances, which require sufficient depth for a correct dispersion
relationship. This paper presents a new method to tackle this dilemma through
an efficient hybridization of volumetric and surface-based advection-
projection discretizations. We design a hybrid time-stepping algorithm that
combines a FLIP domain and an adaptively remeshed Boundary Element Method
(BEM) domain for the incompressible Euler equations. The resulting framework
captures the detailed water splashes near moving objects with the FLIP method,
and produces convincing water waves with correct dispersion relationships at
modest additional costs.

Boundary Element Method (BEM), Fluid Dynamics, Fluid-Implicit-Particle (FLIP)
Method, Hybrid Discretization, Multiscale Modeling.

††submissionid: 176††copyright: rightsretained††journal: TOG††journalyear:
2021††journalvolume: 40††journalnumber: 6††article: 203††publicationmonth:
12††doi: 10.1145/3478513.3480495††ccs: Computing methodologies Physical
simulation Figure 1. Photorealistic rendering of a simulation comprising three
speedboats riding in a circle on the ocean computed with our hybrid approach
combining FLIP and BEM.

##  1\. Introduction

As the Hong Kong American martial artist and actor Bruce Lee commented, “water
can flow or creep or drip or crash.” From calm waves on a boundless sea to
roaring splashes around a speedboat, intriguing large-scale and small-scale
aspects of dynamic liquids have attracted great attention in computational
physics and computer graphics.

Simulating the high-frequency water splashes near the moving structure is
relatively easy by a volumetric Navier-Stokes solver. One of such solvers is
the popular Fluid-Implicit-Particle (FLIP) method, which was proposed by
Brackbill and Ruppel (1986) as an extension to the Particle-In-Cell (PIC)
method (Harlow, 1964), and introduced to computer graphics by Zhu and Bridson
(2005). The difficulty is to efficiently propagate the water waves away from
the structure. Simulating the wide and deep water with high resolution FLIP
for the low frequency waves seems uneconomic. However, a shallow but wide
simulation compromises quality. The depth of the simulated ocean plays a
critical role in determining the shape of the surface, as pointed out by
Nielsen and Bridson (2011). Their approach recovers what was lost in shallow
simulations using a pre-computed deep coarse simulation. Along similar lines,
adaptive discretizations based on advanced grid structures (Setaluri et al.,
2014; Ando and Batty, 2020; Ibayashi et al., 2020; Ando et al., 2013) have
been extensively explored as an option to alleviate the difficulty in
computational cost and memory brought by the demanding resolution in
volumetric schemes. However, these adaptive spatial discretizations are known
for large dissipation in coarse regions. It is not clear if the long-lasting
wide-spreading water waves can be preserved.

Deviating from the adaptivity-related strategies, reduced fluid models
represented by height field approaches (Tessendorf, 2001; Jeschke et al.,
2018; Jeschke and Wojtan, 2017; Yuksel et al., 2007; Canabal et al., 2016)
have also been widely used in graphics applications for capturing large scale
water waves with correct dispersion relationship. Despite their pronounced
efficiency for wave phenomena, these methods cannot capture the water splashes
near obstacles. It is also non-trivial to two-way couple such approaches with
volumetric Navier-Stokes solvers, due to the difficulty in translating back
and forth between the height-field information and the velocity-pressure
values. Although shallow water equation (SWE) solvers have been coupled with
volumetric solvers (Chentanez et al., 2015; Thuerey et al., 2006) to produce
water splashes, the incorrect dispersion relationship of water waves in deep
water can cause unnatural artifacts.

Taking a different route, we hybridize a local high-resolution volumetric FLIP
discretization and a surface-based BEM ocean discretization for the large,
deep ocean. Through a consistent advection-projection discretization of the
incompressible Euler equations, we develop a natural way to synchronize
velocity-pressure information between the two domains. With a surface
representation of water body, we avoid storing any interior fluid quantities
in a large and deep water body, while the surface mesh can be further adapted
by criteria such as wave shape and velocity. For water splashes, we adopt a
massively parallel FLIP simulation in a small region near the structure and
ocean surface to enable intricate close-shot water dynamics.

Our method is inspired by the following observations. First, the
investigations by Da et al. (2016a) demonstrated the potential of BEM for
reproducing complex fluid features with degrees of freedoms only on a surface
geometry. Although its demos are at small scale, we speculate that it can
produce correct wave dynamics (and it does as shown in the result section) in
large scale simulation, and be most efficient compared to volumetric methods
for liquid with a large volume/surface ratio such as a large ocean tank.
Second, the traditionally challenging aspects of BEM such as solid boundary-
handling can be effectively resolved by hybridizing it with a general-purpose
volumetric discretization (the FLIP method), making the whole system robust.
Finally, the fact that both the BEM and the FLIP method adopt 3D velocity
representations of fluid quantities makes the information transfer between
them very natural.

### Specific Contributions

To summarize, we design a new hybrid scheme to efficiently extend water waves
generated in a high resolution FLIP simulation. Our method consists of FLIP to
BEM coupling where BEM receives correction of volumetric simulation, and BEM
to FLIP coupling, where BEM provides the ambient flux of water for the FLIP
simulation, similar to the fluxed animated boundary method (Stomakhin and
Selle, 2017). Our method captures both splashes near structures in multiple
moving FLIP domains and long spreading waves on a vast and deep ocean. The
water wakes extended by our hybrid approach match oceanography studies (Rabaud
and Moisy, 2013), and brute-force FLIP simulation up to 1.7 billion particles.
It also has performance advantages compared to the narrow-band FLIP method
(Ferstl et al., 2016) and the adaptive octree solver (Goldade et al., 2019)
implemented in the commercial software Houdini (Side Effects Software, 2021).

##  2\. Related Work

#### Boundary Element Methods in Graphics

The boundary element method (BEM) has a wide range of applications in computer
graphics. James and Pai (1999) simulated deformable objects interactively with
the BEM for linear elasticity. More complex solid behaviors like fracture can
also be captured and accelerated through the BEM (Hahn and Wojtan, 2016). The
BEM is adopted for simulating various fluid phenomena including large scale
deep ocean waves (Keeler and Bridson, 2015), fluid chains, water droplets (Da
et al., 2016a), and even magnetic fluids (Huang and Michels, 2020). Besides,
BEM has also been adopted in sound rendering (Zheng and James, 2010) and
geometry processing (Solomon et al., 2017).

#### Hybrid Lagrangian-Eulerian Simulation

Hybrid fluid solvers combining Lagrangian and Eulerian perspectives have
gained popularity in recent years. The Particle-In-Cell (PIC) method (Harlow
and Welch, 1965), as a pioneering work of this kind, is limited by its
excessive numerical dissipation. The FLIP method (Brackbill and Ruppel, 1986;
Zhu and Bridson, 2005) effectively improves this issue by only interpolating
grid changes back to particles. More recent approaches such as the APIC method
(Jiang et al., 2015) and the PolyPIC method (Fu et al., 2017) were proposed to
reduce noise and instability of the FLIP method while keeping numerical
dissipation during particle grid transfer under control. By coupling an
Eulerian fluid solver with a narrow band of FLIP particles, Ferstl et al.
(2016) greatly reduced the particle counts while maintaining an
indistinguishable result from a pure FLIP simulation. Losasso et al. (2008)
incorporated Smoothed Particle Hydrodynamics (SPH) within a grid-based
simulation. Domain decomposition approaches (Golas et al., 2012) effectively
reduce the memory and computational costs by combining Eulerian simulations
with vortex particles. Recently, Yue et al. (2018) coupled material point and
discrete element methods to simulate granular material.

#### 2D Methods

SWE methods have been adopted due to their simplicity. Thuerey et al. (2007)
enhanced this model with the capability of handling overturning waves.
Solenthaler et al. (2011) used SPH particles to the SWE. Procedural waves
(Tessendorf, 2001) are effective on simulating vast oceans, however it is non-
trivial to handle boundaries. Procedural waves could be enhanced by satisfying
solid boundary conditions (Jeschke et al., 2020), handling moving obstacles
with wavelet discretization (Jeschke et al., 2018), handling static obstacles
(Jeschke and Wojtan, 2015), or adding high frequency details while matching
the given coarse waves (Nielsen et al., 2013). Schreck et al. (2019) simulated
water waves through fundamental solutions which do not require grids or
meshes. However, in their boat wake scene, the wake does not show the feather-
like pattern in reality. Lagrangian wave particles simulate water waves with
propagating particles (Yuksel et al., 2007; Jeschke and Wojtan, 2017). Canabal
et al. (2016) applied pyramids of convolution kernels with finite support to
2D height fields to get realistic wave dispersion relationships. The amplitude
of the wave animation is given as a parameter tuned by artists. 2D methods are
limited to only capture waves.

#### 2D-3D Coupled Methods

To capture water splashes, coupling 2D methods with 3D methods is a natural
choice. Chentanez et al. (2015) coupled particles, Eulerian grid solver, and
SWE. Thuerey et al. (2006) coupled the SWE with a 3D solver employing a
Lattice Boltzmann Method (LBM). However, as commented by Chentanez et al. “the
SWE and 3D grids simulate inherently different physical phenomena”. The
waveform and propagation speed is visibly different from that of a 3D solver.

#### Spatial Adaptivity

Fluid simulations in large domains such as oceans could be troublesome with
uniform volumetric solvers due to the high computational costs and limited
memory capacities. Adaptive methods such as octrees (Losasso et al., 2004;
Ando and Batty, 2020; Hong, 2009; Aanjaneya et al., 2017) are usually
dissipative in coarse regions. It is not sure if the water waves could be
preserved away from the fine region. Tall cell methods (Irving et al., 2006;
Chentanez and Müller, 2011) were initially designed for shallow water regime.
As reported by Irving et al. (2006), the performance gain due to use of tall
cells is only about a factor of two compared to uniform methods. It offers no
adaptivity along horizontal directions, so simulating wide water bodies is
still challenging for tall-cells methods. Chimera grid methods (English et
al., 2013) require sophisticated mesh handling, and suffer from inefficient
pressure projection. The far-field grid method (Zhu et al., 2013) is suitable
for propagating waves from a single domain, however, supporting multiple
moving domains is difficult. Nielsen and Bridson (2016) adopted adaptive tile
tree for FLIP simulation. Goldade et al. (2019) incorporated an adaptive
octree for viscosity and pressure projection.

#### Guided Simulation

Nielsen and Bridson (2011) proposed guided shapes for high resolution
simulations to reduce the design cycles. Later, Nielsen et al. (2017) further
improved this method by localizing the guided shapes only requiring the
surface mesh. Bojsen-Hansen and Wojtan (2016) proposed a way for achieving
generalized non-reflecting boundary conditions for localized re-simulation.
Absorbing boundary conditions has also been investigated (Söderström et al.,
2010). Similar to our approach, flux boundary method (Stomakhin and Selle,
2017) coupled outer prescribed ocean motion with FLIP/MPM solver through
boundary flux, while it is non-trivial to achieve two way coupling with their
method.

#### Lagrangian Methods

Lagrangian particles such as SPH are also adopted for fluid simulation (Ihmsen
et al., 2014; Huang et al., 2019). Secondary effects such bubble, foam, and
sprays (Ihmsen et al., 2012) can visually enhance the look of liquid in post
processing. Recently, adaptivity has also been investigated in the context of
SPH particles (Winchenbach et al., 2017).

#### Engineering Paradigms

Early work from Longuet-Higgins and Cokelet (1976) utilized the BEM to
numerically compute steep water waves. Higher order time-stepping techniques
with BEM are essential for resolving more general wave phenomena (Grilli et
al., 1989). The BEM could also be combined with the Finite Element Method
(FEM) and be applied for analyzing the nonlinear interactions between waves
and bodies (Wu and Taylor, 2003). More applications and experiments like
propeller blades (Kinnas and Hsin, 1992), added resistance of boats in waves
(Chen et al., 2018), and bubble dynamics (Zhang and Liu, 2015) have been
effectively conducted using the BEM.

##  3\. Method

Figure 2. A moving cage follows the movement of a boat on the ocean. The large
body of water is solved approximately with the BEM (only top layer rendered),
and the resulting water flow enables high resolution simulation with the FLIP
method inside the cage.

The core idea can be intuitively explained by the following imaginary scene: a
boat is moving on the ocean, while there is a moving cage containing the boat,
water, and air as illustrated in Figure 2. We use a large closed mesh to
represent the water of the ocean. Initially there were six faces. The top face
is free surface while other five faces are kinematic solid boundaries. The six
faces are further tessellated into triangles for the use of BEM simulation (Da
et al., 2016a). In Figure 2, we only render the top triangles to avoid visual
distraction. The FLIP method is used exclusively inside the cage. Liquid-boat
interactions are handled by the FLIP method, and the BEM vertices inside the
cage follow the motion of the FLIP liquid surface. The BEM in return provides
the water flux below free surface on the cage boundary for the FLIP
simulation. FLIP influences the BEM only by changing its vertex velocity, and
the BEM influences FLIP only by providing the flux boundary condition below
water. It is similar to the FAB method (Stomakhin and Selle, 2017), except the
flux boundary condition in the FAB method is prescribed by animators, while
the flux boundary condition here is provided by the BEM, which is not
prescribed but driven by the FLIP simulation. This allows waves generated by
solids to propagate out of the FLIP domain, and even to another FLIP domain
via the BEM mesh.

We briefly cover both FLIP and BEM, and then explain our main contribution,
the coupling workflow between them. The entire workflow is summarized in
Algorithm 1 and illustrated in Figure 3.

###  3.1. FLIP and BEM

The governing equation system of our work is the incompressible Euler equation
(Bridson and Müller-Fischer, 2007):

(1) $\displaystyle\frac{D\bm{u}}{Dt}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho}\nabla P+\bm{g}\,,$ (2)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=$ $\displaystyle 0\,,$

where $D/Dt:=\partial/\partial t+\bm{u}\cdot\nabla$ is the material
derivative, $\bm{u}$ is the velocity field, $\rho$ is the constant liquid
density, $P$ is the pressure, and $\bm{g}$ is the gravitational acceleration.

In both, BEM and FLIP, the equation is solved by operator splitting in an
advection-projection style. Assume, in the beginning of the time step, the
velocity is $\bm{u}^{-}=\bm{u}(t)$. After the advection, the velocity becomes
$\bm{u}$. Finally, after the projection, the velocity becomes
$\bm{u}^{+}=\bm{u}(t+\Delta t)$, where $\Delta t$ is the time step size.

In the advection step, we solve the advection equation

(3) $\frac{D\bm{u}}{Dt}=0$

to obtain $\bm{u}$ from $\bm{u}^{-}$, and then add gravity to it. In the
pressure projection step, the incompressibility constraint
$\nabla\cdot\bm{u}=0$ is enforced by solving ${\partial\bm{u}}/{\partial
t}=-{\rho}^{-1}\nabla P$. The associated Poisson’s equation is

$\displaystyle\nabla^{2}P$ $\displaystyle=$ $\displaystyle\frac{\rho}{\Delta
t}\nabla\cdot\bm{u}\,,$ (4) $\displaystyle P|_{\Gamma_{\text{D}}}$
$\displaystyle=$ $\displaystyle P_{\text{bc}}\,,$
$\displaystyle\frac{dP}{dn}|_{\Gamma_{\text{N}}}$ $\displaystyle=$
$\displaystyle\frac{\rho}{\Delta
t}(\bm{u}-\bm{u}_{\text{solid}})\cdot\bm{n}\,,$

where $\Gamma_{\text{D}}$ and $\Gamma_{\text{N}}$ are liquid-air (Dirichlet)
and liquid-solid (Neumann) boundaries respectively. $P_{\text{bc}}$ is the
Dirichlet pressure boundary condition which usually consists of surface
tension, $\bm{u}_{\text{solid}}$ is the desired liquid velocity on the solid
boundary, which can introduce inflow and outflow, and $\bm{n}$ is the external
normal vector of the liquid domain. Eq. (4) is solved to get the pressure $P$,
and the velocity is then updated by

(5) $\bm{u}(t+\Delta t)=\bm{u}^{+}=\bm{u}-\frac{\Delta t}{\rho}\nabla P\,.$

The pressure projection Eq. (4) is solved differently in the BEM and FLIP. For
the FLIP simulation, we generally follow Zhu and Bridson (2005). In the
pressure projection step, we use the approach of Ng et al. (2009) to handle
the liquid-solid boundary, and refer to Gibou et al. (2002) for the free-
surface boundary condition. Since we focus on large scenes, we neglect the
surface tension.

In short, the BEM simulates the motion of liquid by evolving triangle meshes.
Each vertex carries position and velocity information of the liquid surface.
In the advection step, each vertex moves with corresponding velocity to the
new position. However, when these vertices move, collisions between points to
triangles, and edges to edges happen. Furthermore, when two droplets collide,
topology changes are required to join two disjoint droplet meshes together.
Such collision handling and topology changes are handled by a mesh-based
surface tracking program LosTopos (Da et al., 2014). The input of LosTopos is
the original vertex positions, triangle connectivity, and the target vertex
positions. It resolves collisions, topology changes, and outputs a new set of
vertex positions and triangle connectivity. The output mesh is intersection-
free and watertight. After the advection, the surface velocity field may not
represent the boundary value of a divergence-free velocity field. Therefore,
it is followed by a projection step consisting of Helmholtz decomposition
(H.D.) and solving a Laplace equation with mixed boundary condition. We refer
the reader to the work of Da et al. (2016a), and the work of Huang and Michels
(2020) for details. We prefer the H.D. in the implementation of Da et al.
(2016a) because the look is closer to our reference 3D simulation, see Figure
9(a) and 9(b). The theoretical analysis of such preference is presented in the
supplemental material.

Input: FLIP particles and BEM mesh at time $t$.

Output: FLIP particles and BEM mesh at time $t+\Delta t_{f}$.
$t_{\text{BEM}}\leftarrow 0$;

while  _$t_{\text{BEM}} <\Delta t_{\text{frame}}$_ do

$\Delta t_{\text{BEM}}\leftarrow$ min(_$\Delta x_{\text{BEM}}$ /
max(_$v_{\text{BEM}}$_), $\Delta t_{\text{frame}}-t_{\text{BEM}}$_);

$t_{\text{BEM}}\leftarrow t_{\text{BEM}}+\Delta t_{\text{BEM}}$;

/* Figure 3, Plot (2), Figure 4 (left), Section 3.2 */

FLIPtoBEMCoupling(_$\Delta t_{\text{BEM}}$ , FLIP SDF and velocity_);

BEMAdvection(_$\Delta t_{\text{BEM}}$_);

/* Figure 3, Plot (3) */

BEMProjection(_$\Delta t_{\text{BEM}}$_);

end while

$t_{\text{FLIP}}\leftarrow 0$;

while  _$t_{\text{FLIP}} <\Delta t_{\text{frame}}$_ do

$\Delta t_{\text{CFL}}=$ $\text{cfl}()\Delta x_{\text{FLIP}}$ /
max(_$v_{\text{FLIP}}$_);

$\Delta t_{\text{FLIP}}\leftarrow$ min(_$\Delta t_{\text{CFL}}$ , $\Delta
t_{\text{frame}}-t_{\text{FLIP}}$_);

$t_{\text{FLIP}}\leftarrow t_{\text{FLIP}}+\Delta t_{\text{FLIP}}$;

/* Figure 3, Plot (4) */

FLIPAdvection(_$\Delta t_{\text{FLIP}}$_);

/* Figure 4 (right), Section 3.3 */

BEMtoFLIPCoupling(_BEM mesh and velocity_);

/* Figure 3, Plot (5) */

FLIPProjection(_$\Delta t_{\text{FLIP}}$_);

end while

ALGORITHM 1 A single frame is simulated using our coupled BEM and FLIP
approach. A time step size of $\Delta t_{\text{frame}}$ is used, $\Delta
x_{\text{BEM}}$ denotes the minimum edge length of the BEM mesh,
$v_{\text{BEM}}$ indicates the BEM vertex velocity, $\Delta x_{\text{FLIP}}$
denotes the FLIP voxel size, and the FLIP velocity is given by
$v_{\text{FLIP}}$. The CFL number is given by cfl(). Extra coupling steps are
written in blue.
\begin{overpic}[width=424.94574pt]{graphics_final/whole_flow_chart_notext.pdf}
\put(0.0,30.0){(1)} \put(30.0,30.0){BEM} \put(34.5,30.0){(2)}
\put(42.0,30.0){Guided Advection} \put(69.0,30.0){(3)}
\put(80.0,30.0){Projection} \put(30.0,14.0){FLIP} \put(34.5,14.0){(4)}
\put(45.0,14.0){Advection} \put(69.0,14.0){(5)} \put(76.0,14.0){Projection
with Flow} \put(7.0,3.0){Next frame FLIP boundary} \put(7.0,4.73){FLIP
boundary} \put(7.0,6.45){BEM mesh} \put(7.0,8.17){Interpolated BEM velocity}
\put(7.0,9.89){FLIP velocity} \put(7.0,11.61){BEM velocity} \end{overpic}
Figure 3. Flow chart of a single time step. (1) In this scene a water ball
impacts the surface causing it to bend inward. The light blue region is the
FLIP liquid domain. The solid dark blue and red lines are the boundary of the
FLIP liquid and the BEM surface mesh respectively. The dashed blue line
indicates the FLIP free surface after advection. The dark blue and red arrows
indicate the liquid velocity on the boundary of the FLIP and BEM simulation
respectively. The black arrows are the FLIP boundary velocity interpolated
from the BEM mesh. (2) BEM vertices move to its new positions guided by the
FLIP free surface and velocity; see Section 3.2 and Figure 4 (left). (3) After
the BEM projection step, the velocity field becomes divergence free again. (4)
FLIP advection. (5) The boundary velocity of the FLIP domain is interpolated
from the BEM surface velocity field. It is then plugged in the projection step
of FLIP to make the velocity field divergence free; see Section 3.3. This
concludes a single time step.

###  3.2. FLIP Influence on the BEM

\begin{overpic}[width=203.80193pt]{graphics_final/new_BEM_correction_illustration.pdf}
\put(10.0,47.0){BEM vertex} \put(9.0,41.0){$x_{0}$} \put(8.0,23.0){$x_{1}$}
\put(22.0,6.0){$x_{2}$} \put(22.0,1.0){$v_{\text{FLIP}}$}
\put(36.0,41.0){$x_{3}$} \put(21.0,15.0){$x_{4}$} \put(23.0,31.0){$x_{5}$}
\put(22.0,23.0){$x_{\text{final}}$} \put(2.0,14.0){$x_{\text{sample}}$}
\put(30.0,10.0){FLIP surface} \put(58.0,5.0){Flow} \put(75.0,5.0){BEM domain}
\put(68.0,11.0){{\color[rgb]{0.21484375,0.44140625,0.78515625}\definecolor[named]{pgfstrokecolor}{rgb}{0.21484375,0.44140625,0.78515625}Fill
zone}} \put(70.0,28.0){Boat} \put(65.0,34.0){FLIP domain}
\put(63.0,41.0){{\color[rgb]{1,0.66796875,0.66796875}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.66796875,0.66796875}Sink
zone}} \end{overpic} Figure 4. Illustration of the advection step in the
coupling mechanism between BEM and FLIP. The details are explained in Section
3.2 and 3.3 . Left: A BEM vertex at $x_{0}$ moves under the guidance of FLIP
to its final destination at $x_{\text{final}}$. Right: A boat is moving on the
ocean represented by the BEM domain. A buffer region around the FILP domain is
divided into the sink zone where particles are removed, and the fill zone
where new particles are generated.

Figure 5. Top: The BEM mesh moves under the guidance of FLIP particles,
closely matching its surface. The incoming boat is removed for better
visualization of the bottom FLIP surface. With reference coordinate system
fixed at the boat, white color indicates slow velocity while blue indicates
high velocity. Bottom: Particles move to the left, exiting the FLIP domain.
Particles in the buffer zone and above the BEM mesh are removed, while
particles below the BEM mesh are filled to keep constant density.

This section corresponds to the guided advection step in Figure 3. The BEM
takes information from the FLIP simulation by tracking the water flow on the
liquid-air and liquid-solid interfaces of the FLIP system. However, just
advecting BEM vertices with FLIP velocity field may leads to a drift of free
surface. Therefore we must consider both the flow velocity and liquid
interface.

We expect the FLIP to BEM coupling to be non-intrusive to the existing BEM
program, so we decide to achieve the goal by modifying the velocity of each
vertex before the BEM time stepping function, instead of directly matching the
interface position and flow velocity separately inside the advection function
of the BEM.

There are splashes in the FLIP simulation. Letting the BEM track all these
detailed droplets is unnecessary because those droplets will eventually merge
into the main body water. When they merge, the BEM is going to track the main
surface again. Therefore, we let the BEM track a smoothed FLIP liquid signed
distance function (SDF).

The detailed algorithm is as follows. The key steps are illustrated in Figure
4. Snapshots from actual simulations are shown in Figure 5. We use the
smoothed SDF of the FLIP liquid and velocity field at the beginning of current
frame as a guidance to find the target position that a BEM vertex is supposed
to be at the end of a substep. Based on that we modify the velocity.

First, given the current BEM vertex position $x_{0}$, we search from a point
below $x_{0}$ and inside the FLIP fluid, stepping $0.5\Delta x_{\text{FLIP}}$
upward until the sampled liquid SDF becomes positive, then based on the
sampled positive SDF we move downward to get the sub-voxel position of the
interface $x_{1}$. If there is no intersection $x_{1}$, then all following
operations are skipped.

In the current cut-cell formulation of Neumann boundary conditions in the FLIP
(Ng et al., 2009), we find it helpful to produce water splashes in front of
the boat if we mix the liquid velocity with solid velocity proportional to the
solid fraction of the corresponding voxel face. For example, if the voxel face
of a velocity component has $30\%$ of solid, with an artificial friction
coefficient of 0.8, and the boat is static, the blended component velocity
would be $76\%$ of the original value. However, the side effect is that BEM
takes such low velocity as a sign of large resistance, producing exaggerated
wake amplitudes (Figure 10). To alleviate this side effect, we sample the FLIP
velocity not directly at the surface point $x_{1}$, but a few voxels below:

(6) $\displaystyle x_{\text{sample}}$ $\displaystyle=$ $\displaystyle
x_{1}-(0,\alpha\Delta x_{\text{FLIP}},0)\,,$ (7) $\displaystyle
v_{\text{FLIP}}$ $\displaystyle=$
$\displaystyle\text{VelocitySample}(\text{FLIP},x_{\text{sample}})\,.$

The BEM still tracks the correct free surface, but feels less resistance due
to faster flow velocity sampled from the FLIP. The effects of sampling depth
$\alpha$ will be investigated in Section 5.5 and in Figure 10.

The sampled velocity $v_{\text{FLIP}}$ determines the ideal target location of
the BEM vertex $x_{2}$ at the end of the frame:

(8) $x_{2}=x_{1}+\Delta t_{\text{frame}}v_{\text{FLIP}}\,.$

A problem with setting $x_{2}$ as the substep target position for the BEM
vertex is that it would attempt to have a high velocity to travel very long
distance ($x_{2}-x_{0}$) in a short period of time if there is initial surface
mismatch or substep time $\Delta t_{\text{BEM}}$ is smaller than frame time
$\Delta t_{\text{frame}}$. Therefore, we limit the distance that a BEM vertex
can travel by its original velocity $v_{\text{BEM}}$ and substep time $\Delta
t_{\text{BEM}}$ to get the limited surface target $x_{4}$:

(9) $\displaystyle x_{3}$ $\displaystyle=$ $\displaystyle x_{0}+\Delta
t_{\text{BEM}}v_{\text{BEM}}\,,$ (10) $\displaystyle x_{4}$ $\displaystyle=$
$\displaystyle
x_{0}+\frac{x_{2}-x_{0}}{|x_{2}-x_{0}|}\text{min}(|x_{2}-x_{0}|,|x_{3}-x_{0}|)\,.$

However, if all BEM vertices are initially static, they would not move at all.
Therefore we need to blend the velocity with the FLIP flow velocity
$v_{\text{FLIP}}$. The target position with pure FLIP flow velocity is
$x_{5}$:

(11) $x_{5}=x_{0}+\Delta t_{\text{BEM}}v_{\text{FLIP}}\,.$

The final substep target position of the BEM vertex is a linear mixture of
free surface target $x_{4}$ and FLIP flow target $x_{5}$. We just take their
mid point:

(12) $\displaystyle x_{\text{final}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(x_{4}+x_{5})\,.$

Finally, the vertex velocity is obtained:

(13) $v_{\text{BEM}}=\frac{x_{\text{final}}-x_{0}}{\Delta t_{\text{BEM}}}\,.$

Although the SDF and velocity of the FLIP are taken at the beginning of the
frame, the surface target position $x_{2}$ is at the end of the frame as an
approximation of where the FLIP surface is supposed to be. An alternative way
is to advect the SDF and velocity to the end of BEM substep time, and design
another workflow, but it is more costly. When it is not the first substep, the
projection from $x_{0}$ to $x_{1}$ actually do not happen at the beginning of
the frame. To force it to happen at the frame beginning, one can estimate
where the vertex would be by backtracking with the vertex velocity backward to
the frame beginning, find the vertical projection point $x_{1}$, and estimate
the target position $x_{2}$ at the end of the frame. However, we did not
observe visual differences in this more intricate approach. Also most of the
time there is only one BEM substep.

###  3.3. BEM Influence on the FLIP

The BEM provides the water flow boundary condition of the FLIP domain. The
difficulty is to know the velocity inside the BEM mesh. If this interpolation
problem is to be solved by potential flow or Stokes flow (Bhattacharya et al.,
2012), the whole BEM domain would be voxelized and equations are inefficiently
solved. Finding the nearest point velocity on the surface seems to be the
simple solution, but the velocity field is not physical, see Figure 8.
Therefore we propose to use boundary integral over BEM triangles to get the
velocity on the FLIP Neumann boundary inside BEM mesh:

(14)
$\bm{u}(x)=\nabla\int_{\Gamma}\frac{\bm{n}(y)\cdot\bm{u}(y)}{4\pi|x-y|}ds_{y}-\nabla\times\int_{\Gamma}\frac{\bm{n}(y)\times\bm{u}(y)}{4\pi|x-y|}ds_{y}\,,$

where $x$ is inside the BEM domain, and $y$ is on the BEM boundary $\Gamma$.
If we introduce two auxiliary variables $c$ and $\bm{j}$ (similar to charges
and currents in electrodynamics), the above equation becomes:

(15) $\displaystyle c(y)$ $\displaystyle=$
$\displaystyle\bm{n}(y)\cdot\bm{u}(y)\,,$ (16) $\displaystyle\bm{j}(y)$
$\displaystyle=$ $\displaystyle\bm{n}(y)\times\bm{u}(y)\,,$ (17)
$\displaystyle\bm{u}(x)$ $\displaystyle=$
$\displaystyle\nabla\int_{\Gamma}\frac{c(y)}{4\pi|x-y|}ds_{y}-\nabla\times\int_{\Gamma}\frac{\bm{j}(y)}{4\pi|x-y|}ds_{y}\,,$
(18) $\displaystyle=$
$\displaystyle\nabla\int_{\Gamma}\frac{c(y)}{4\pi|x-y|}ds_{y}+\int_{\Gamma}\bm{j}(y)\times\nabla\frac{1}{4\pi|x-y|}ds_{y}\,.$

Since $x$ is not on the boundary $\Gamma$, $|x-y|\neq 0$. The above integral
is not singular. The integral is calculated as the summation of integrals of
each triangle on the BEM mesh. On each triangle, the normal direction $\bm{n}$
is constant and the velocity $\bm{u}$ is the barycentric interpolation of
vertex values. Accordingly $c(y)$ and $\bm{j}(y)$ on the triangle are
calculated based on the interpolated velocity. On each triangle, we select
three symmetric quadrature points, whose barycentric coordinates are
$2/3,1/6,1/6$ with rotation, and weights are just one third of the triangle
area.

The interpolation point $x$ can still be close to the quadrature point which
introduces error. The easy solution is to use nearest neighbor interpolation,
or one can apply the more expensive analytic integration of the fundamental
solution $1/|x-y|$ over

triangles, as adopted by Huang and Michels (2020). The brute force integration
for each position on the FLIP boundary is expensive. We take the advantage
that the interpolated velocity field inside BEM mesh is smooth. We create a
narrow band velocity volume which typically has grid resolution of $4\Delta
x_{\text{FLIP}}$ around the FLIP boundary. The interpolated velocity based on
boundary integral is zero if $x$ is outside the BEM mesh. Therefore when
creating this narrow band velocity volume, we only sample velocity inside the
BEM mesh, and extrapolate one layer to enhance coverage. Otherwise the
velocity near BEM free surface would be affected by the zero velocity sample
outside. In addition to boundary velocity, we also need to add and remove
particles. Otherwise when the water flows in, FLIP particles left the influx
Neumann boundary, and there is a gap between particles and Neumann boundary.
Therefore we just need to fill the gap. The gap size is related to the CFL
condition, so it is usually a few voxels after advection. Inside the Neumann
boundary, we set a fill zone, where particles are filled. Since our maximum
CFL is 4, we use 4 voxels wide fill zone. A voxel in the fill zone and inside
the BEM mesh are filled with particles, whose initial velocity is given by the
narrow band velocity volume obtained from boundary integral. We also set a
sink zone where particles are removed. The sink zone is similarly defined
except it is outside of the BEM mesh. See Figure 4 and 5.

##  4\. Implementation

The BEM part is based on the open source code of Da et al. (2016b) with CUDA
acceleration provided by Huang and Michels (2020). The FLIP code is based on
the open source implementation of Batty et al. (2007) and Ng et al. (2009)
provided by Batty (2018). However, we re-wrote the whole FLIP simulation with
OpenVDB (Museth, 2013) to support sparse data structures.

During the guided advection step in the BEM, the liquid SDF of FLIP is
represented as an OpenVDB float grid. We first re-initialize the liquid SDF in
FLIP simulation with OpenVDB function “levelSetRebuild”, since only in the
one-voxel neighborhood of a particle we have a valid SDF, while the following
Gaussian filter needs wider bandwidth. We then apply three passes of one-
voxel-wide Gaussian filters provided by OpenVDB to smooth the liquid SDF. As
illustrated in Figure 4 (left), we need to find the BEM vertex vertical
projection point on the FLIP surface. We search from the bottom of the FLIP
domain, and step upward until the sampled liquid SDF becomes positive. The
search step size is half of the FLIP voxel size. The projection step of the
BEM is a Galerkin boundary element method as described and implemented by
Huang and Michels (2020).

In the grid-to-particle step of FLIP, we use 0.95 FLIP velocity mixed with
0.05 PIC velocity. In the particle-to-grid step we transfer not only the
momentum but also the liquid’s SDF. Each FLIP particle is a liquid ball, whose
diameter is 1.01 times half of the voxel diagonal.

For moving FLIP domains, we update the location of the fill and sink zones at
the beginning of each sub-step before particle advection. The sub-step size
must be small enough so no particles moves outside the updated fill zone. All
air gaps are immediately filled by particles after the advection step. In the
fill and remove particles step, we convert the BEM mesh to a SDF with
functions provided by OpenVDB. When we fill new particles, we must ensure that
after the new particles turned into a smoothed SDF, the smoothed FLIP SDF
should coincide with the BEM SDF, so that the water level is not altered by
coupling. We found generating a particle if its center is slightly ($0.7$ of
the particle radius) interior to the BEM surface reproduces the same SDF after
smoothing. We remove all FLIP particles in the sink zone (domain buffer zone
where the BEM SDF is greater than the -0.7 radius threshold). Around the fill
zone (domain buffer zone with BEM mesh and SDF smaller than zero) we create a
narrow band velocity volume and sample the interpolated velocity every four
voxels according to Section 3.3. We count the particle numbers in each voxel
in the fill zone. We generate new particles at random positions in a voxel if
the count is less than 8, and remove extra particles if it is larger than 16.
The newly generated particles in the fill zone takes velocities linearly
interpolated from the narrow band velocity volume. This approach saves the
expensive boundary integral in Eq. (14).

The FLIP domain boundary velocity is also linearly interpolated from above
narrow band velocity volume. In the pressure projection part we generally
follow Weber et al. (2015), except we use unsmooth aggregation to build the
multigrid hierarchy algebraically and replace their V cycle by the W cycle. We
apply three Red-Black Gauss-Seidel iterations before restriction and after
prolongation in the W cycle. To define the restriction operator, each eight
neighbor voxels at current level are averaged to become a larger voxel at the
coarser level. The coarser level matrix is scaled by 0.5 according to Stüben
(2001) (see Section 6.4 therein).

##  5\. Results

Table 1. Overview of the different scenes presented in this paper. The resolution $\Delta x$ is stated in meter, and the time step size $\Delta t$ and per-invocation run time $T$ in seconds. We list the invocation counts in the parenthesis for each type of run time since adaptive sub-steps happen. The coupling is carried out once per frame. Eight particles occupy a single voxel when using FLIP. Machine A has 2 $\times$12-core Intel Xeon<EMAIL_ADDRESS>and NVIDIA RTX 2080Ti. Machine B has 2$\times$64-core AMD EPYC<EMAIL_ADDRESS>Machine C has 2 $\times$16-core Intel Gold<EMAIL_ADDRESS>and RTX2080Ti. Machine D has 2$\times$24-core Intel Platinum<EMAIL_ADDRESS>and NVIDIA V100. Scene | Figure | $N_{\text{Tri}}$ | $N_{\text{Particle}}$ | $\Delta x_{\text{BEM}}$ | $\Delta x_{\text{FLIP}}$ | $\Delta t$ | $T_{\text{Frame}}$ | $T_{\text{BEM}}$ | $T_{\text{FLIP}}$ | $T_{\text{Coupling}}$ | Machine
---|---|---|---|---|---|---|---|---|---|---|---
Dam Break | 6 | 20K | / | 0.03 | / | 1/2000 | 4.7 (4000) | 3.8 (4925) | / | / | A
Crown Splash | 9(a) | 24K | / | 0.02 | / | 1/60 | 6.8 (600) | 5.2 (786) | / | / | A
9(b) | 33K | / | 0.02 | / | 1/60 | 8.2 (600) | 5.3 (932) | / | / | A
9(c) | 20K | 5M | 0.02 | 0.01 | 1/60 | 6.3 (600) | 3.7 (717) | 0.7 (999) | 0.6 | A
9(d) | / | 5M | / | 0.01 | 1/60 | 2.0 (600) | / | 0.7 (1707) | / | A
9(e) | / | 15M | / | 0.01 | 1/60 | 4.3 (600) | / | 1.62 (1598) | / | A
Boat (Straight, |  10f | / | 1705M | / | 0.14 | 1/24 | 70.5 (480) | / | 70.5 (480) | / | B
Half Resolution) |  10g | / | 320M | / | 0.14 | 1/24 | 21.6 (480) | / | 21.6 (480) | / | B
Boat (Straight) |  10d | 65K | 99M | 0.7 | 0.07 | 1/24 | 49.3 (480) | 10.1 (480) | 10.6 (1415) | 7.1 | C
10h | 55K | 99M | 0.7 | 0.07 | 1/24 | 42.4 (480) | 7.4 (480) | 9.3 (1414) | 6.8 | A
15a | 111K | 21M | 0.14 | 0.07 | 1/100 | 59.7 (500) | 55.3 (501) | 2.4 (500) | 1.8 | A
15b | / | 210M | / | 0.07 | 1/100 | 18.2 (500) | / | 18.2 (500) | / | A
Boat (Wavy) | 16 | 89K | 230M | 0.7 | 0.07 | 1/24 | 105.1 (480) | 15.8 (923) | 22.5 (1404) | 9.8 | C
Boats (Cycling) |  1, 17 | 170K | 634M | 0.7 | 0.07 | 1/24 | 248.9 (600) | 41.2 (600) | 64.3 (1766) | 15.8 | D

In this section, we present examples to evaluate our hybrid method. In all
experiments, we use eight particles per voxel in the FLIP simulation, and
reduce the maximum relative divergence error in the pressure equation to
$10^{-7}$. First a BEM dam break example shows that the adopted BEM
formulation captures 3D fluid physics. Then we use a reference FLIP simulation
to drive both the shallow water equation (Chentanez et al., 2015) and BEM, to
show that the BEM can extend water waves with the correct shape and speed
while SWE cannot. After we test the FLIP to BEM coupling, we show our BEM to
FLIP coupling based on the boundary integral is compatible with Stokes’s wave
theory (Stomakhin and Selle, 2017). Next we validate the two-way coupling by
extending a small FLIP region with BEM to a larger tank, whose motion
qualitatively matches the reference FLIP solution in the larger tank.
Furthermore, we evaluate our method where FLIP handles liquid-solid collisions
to complement the weakness of BEM. We explore the choice of parameters on the
shape of the boat wake pattern, and show the wake pattern can match real world
cases measured on airborne images (Rabaud and Moisy, 2013). Finally we show
that our hybrid method supports multiple FLIP domains, and supports secondary
white water simulation due to the horizontal velocity field provided by the
BEM mesh. Relevant parameters and run times are summarized in Table 1.

###  5.1. Dam Break

Figure 6. Dam break simulated by BEM following the setup of Colagrossi and
Landrini (2003). Our choice for the BEM formulation can capture the defining
turning waves of this experiment as well as the turbulent motions of the
surface.

Our first experiment showcases the ability of the surface-only liquids
formulation to capture 3D liquid physics. This distinguishes our BEM solver
from height-field- or wave-equation-based options. We position a
$2\,\textrm{m}\times 1\,\textrm{m}\times 0.2\,\textrm{m}$ fluid volume at the
corner of a $5.366\,\textrm{m}\times 3\,\textrm{m}\times 0.2\,\textrm{m}$
container according to the setup of Colagrossi and Landrini (2003) and let it
flow freely. Our choice of the BEM solver is able to qualitatively match the
shape of their simulation as shown in Figure 6.

###  5.2. FLIP Influencing SWE and BEM

Previous methods on 2D-3D coupling extend the 3D fluid simulation with shallow
water equation (Thuerey et al., 2006; Chentanez et al., 2015). However,
shallow water equation cannot propagate the water wave with the correct shape
and speed. We drop a water ball with 0.7 m radius from 3 m above the ground at
the center of a $20\,\textrm{m}\times 1\,\textrm{m}\times 20\,\textrm{m}$
shallow tank. The SWE, FLIP and BEM simulation cover the same domain, with
$\Delta x=0.1,0.1,0.2\textrm{m}$ respectively. The reference FLIP simulation
drives the SWE simulation within a 6 m-radius circular region at the center
according to (Chentanez et al., 2015), and drives the BEM simulation within a
square region with 6 m edge length at the center. The result is shown in
Figure 7. When waves propagate out of the influence region, the SWE water
waves have incorrect shapes and speeds while BEM waves can qualitatively match
the reference. Note this example already use a shallow tank. The SWE would
perform worse in a deeper tank.

\begin{overpic}[width=203.80193pt]{graphics_final/composed_frame124_with_colorbar.png}
\put(15.0,51.0){SWE Elevation} \put(63.0,51.0){BEM Elevation}
\put(39.0,44.0){FLIP Elevation} \put(10.0,5.0){SWE Error} \put(75.0,5.0){BEM
Error} \put(37.0,10.0){-0.2\,m} \put(56.0,10.0){0.2\,m} \end{overpic} Figure
7. SWE and BEM simulations driven by a reference FLIP simulation. Middle: the
reference FLIP simulation. Top left (right): The SWE (BEM) simulation driven
by the reference FLIP. Blue and red indicate low and high elevation. Bottom
left (right): The same SWE (BEM) simulation color encoded with nearest
distance to the reference FLIP mesh. When waves propagate out of the influence
region, SWE water waves have incorrect shapes and speeds while BEM waves can
qualitatively match the reference.

###  5.3. Boundary Integral Velocity Interpolation

We determine the velocity inside the BEM mesh by computing the boundary
integral of the velocity on the mesh. This approach is compatible with
Stokes’s wave theory adopted in the FAB method (Stomakhin and Selle, 2017). We
start with a mesh that represents the shape of a Stokes wave, with correct
mesh surface velocity field according to Stokes’s theory both on and below the
curvy free surface. The simulation is carried out in a reference frame moving
along with the wave, so the mesh look static over time. We carve out a piece
of mesh and fill with FLIP particles. With the correct flow boundary
condition, the FLIP simulation should reproduce the same Stokes wave. In
Figure 8 we compare the nearest neighbor interpolation, our interpolation, as
well as the velocity directly from Stokes’s theory. Our interpolation method
is able to correctly recover the fluid velocity inside the mesh.

\begin{overpic}[width=203.80193pt]{graphics_final/Nearest_BEM_analytic145.png}
\put(0.0,41.0){Nearest Neighbor} \put(0.0,27.0){Ours} \put(0.0,13.0){Stokes
Theory} \end{overpic} Figure 8. Comparison of velocity interpolation method
for a Stokes wave mesh. Top: Nearest neighbor interpolation. Middle: Our
boundary integral interpolation. Bottom: Reference Stokes wave theory
velocity. The nearest neighbor interpolation does not provide a correct
boundary condition for the FLIP to maintain Stokes’s wave. Our interpolation
method can produce the correct result.

###  5.4. Two-Way Coupling for Domain Extension

(a) Pure BEM (Da et al., 2016a).

(b) Pure Full H.D. BEM (Huang and Michels, 2020).

(c) Our hybrid approach.

(d) Center FLIP without coupling.

(e) Reference FLIP.

Figure 9. Simulation of a water ball dropping into a cylinder pool. The pure
BEM simulation (a) and our hybrid FLIP-BEM simulation (c) with particles and
BEM mesh can match the reference FLIP simulation (e), but the BEM with full
H.D. (b) does not match the reference FLIP simulation (e). The FLIP simulation
without coupling (d) is affected by its boundary.

We validate our two-way coupling scheme in a domain extension problem. A water
ball with 0.1 m radius impacts a $1.2\,\textrm{m}\times 0.4\,\textrm{m}\times
1.2\,\textrm{m}$ square FLIP domain from 1 m above the waterline. The waves
are extended by the BEM, and touch the outer cylinder boundary with 1 m radius
and 0.6 m depth. The reflected waves in BEM propagate back to the small FLIP
domain. Ideally, the motion should match a reference FLIP simulation in the
larger domain. In this example, the minimum edge length of the triangle in the
BEM simulation is $0.02\,\textrm{m}$. The voxel size of the FLIP simulation is
$0.01\,\textrm{m}$.

The FLIP particles are integrated with midpoint rule with CFL condition = 1.
The results are shown in Figure 9.

We prefer the partial H.D. (Helmholtz decomposition). formulation of Da et al.
(2016a) over the full H.D. of Huang and Michels (2020). We show the different
H.D. in Figure 9(a) and 9(b) respectively. The BEM solver with partial H.D. is
able to qualitatively match the dynamic motion of the reference FLIP
simulation. However, the full H.D. leads to quite different results compared
to the reference FLIP simulation.

Since BEM alone already matches the reference FLIP qualitatively, using it to
extend the small FLIP domain can probably produce a matching result as well.
This is indeed the case. As shown in Figure 9(c), the hybrid simulation looks
similar to the reference one. Initially, the BEM mesh is static and the
momentum is transferred from the water drop in the center FLIP simulation. Our
hybrid method allows the water wave to propagate beyond the FLIP box boundary
and to reach the true cylinder boundary. There is visible square pattern
because we render the FLIP particles with their simulation radius, so part of
the sphere intersects the BEM mesh. In Figure 9(c) (second from left), the
interior particles rise slightly higher than the transition zone. This is an
artifact of our weakly coupling scheme. Here, when the water ball impacts the
surface, the velocity interpolated from BEM does not generate enough outflow
of the liquid on the FLIP domain boundary. Therefore, the particles inside the
small FLIP domain was squeezed by the impact and rise up. A strong coupling
scheme may solve this issue, but it may be much more expensive to solve the
system. The BEM vertex inside the FLIP domain is governed by FLIP system with
BEM velocity boundary, while the BEM vertex just outside is governed by BEM
system. The seam indicates how these two systems match. The two system is
different, but most of the time, their difference is small. In contrast, as
shown in Figure 9(d), without the BEM coupling, the FLIP simulation alone is
influenced by the box boundary leading to different fluid motion.

###  5.5. Wake Pattern

\begin{overpic}[width=420.61192pt]{graphics_final/with_friction_01234FLIPrefshallow_nofriction3.png}
\put(0.5,22.0){(a)} \put(13.0,22.0){(b)} \put(25.5,22.0){(c)}
\put(38.0,22.0){(d)} \put(50.5,22.0){(e)} \put(63.0,22.0){(f)}
\put(75.5,22.0){(g)} \put(88.0,22.0){(h)} \end{overpic} Figure 10. Wake
patterns of a sailing boat on a $120\,\textrm{m}\times 240\,\textrm{m}$ domain
with different configurations. Figures a-e correspond to hybrid methods with
$\Delta x_{\text{FLIP}}=0.07\,\textrm{m}$, $\Delta
x_{\text{BEM}}=0.7\,\textrm{m}$, FLIP depth $3.5\,\textrm{m}$, artificial
friction coefficient = 0.8, and BEM samples FLIP velocity 0,1,2,3,4 voxels
below the FLIP surface. Figure f and g correspond to FLIP simulation with $20$
and $3.5\,\textrm{m}$ depth, and $\Delta x_{\text{FLIP}}=0.14\,\textrm{m}$.
Figure h is without artificial friction, and samples three voxels below the
interface. Blue and red colors map to -0.8 m and 0.8 m elevation in all
figures. Black boxes indicate the FLIP regions.

We analyze our hybrid method for simulating the ship/boat wakes. In
particular, we analyze key factors influencing the wake amplitudes, compare
the wake angles to real world measurements, and examine the overturning waves
behind the ship/boat.

The wake amplitudes and wake angle tests are carried out in a
$120\,\textrm{m}\times 20\,\textrm{m}\times 240\,\textrm{m}$ domain. The boat
is about 10 m long. The coupled FLIP domain is $20\,\textrm{m}\times
3.5\,\textrm{m}\times 60\,\textrm{m}$. The water flow velocity is
$10\,\textrm{m/s}.$

The influence of sampling depth $\alpha$ in Eq. (6) is illustrated in Figure
10. In Figure 10a-e, we sample at 0,1,2,3,4 voxels below the FLIP surface
respectively. The artificial friction coefficient is 0.8. In Figure 10f, we
run the reference FLIP simulation, and in Figure 10g, we only simulate a
shallow 3.5 m tank to examine influence of depth on the wake pattern.
Additionally, in Figure 10h we turn off the artificial friction, and sample 3
voxels below. With increasing sampling depth, the amplitudes are closer to
that of the reference FLIP. For the same sampling depth, turning off the
artificial friction also leads to an amplitude closer to the reference as
shown in Figure 10h. The shallow tank simulation in Figure 10g has a different
water wave dispersion relationship compared to Figure 10f.

\begin{overpic}[width=424.94574pt]{graphics_final/friction_vs_nofriction_vs_REF.png}\put(1.0,60.0){(a)}\put(1.0,47.5){(b)}\put(1.0,35.0){(c)}\put(1.0,22.5){(d)}\put(1.0,10.0){(e)}\end{overpic}
Figure 11.  Simulated boat wakes. In realistic renderings, FLIP particles are
converted to mesh, and rendered with BEM mesh. See Section 5.6 for detailed
steps. (a) FLIP particles and BEM mesh from hybrid simulation in Figure 10d.
(b) Realistic rendering of hybrid simulation in Figure 10d. (c) Realistic
rendering of Figure 10h. (d) Realistic rendering of reference FLIP simulation
in Figure 10f. (e) Simulated result with (Canabal et al., 2016). Figure 12.
Real photo of boat wakes at Lyse fjord in Norway, adapted from a photograph by
Edmont, distributed under a CC BY-SA 3.0 license (Edmont, 2009).

In Figure 11, we render some wake patterns to investigate the influence of
wake amplitudes on visual appearance. We choose the hybrid results with 3
voxels of sampling depth, with (Figure 11a,b) and without (Figure 11c)
artificial friction. The reference FLIP simulation is presented in Figure 11d.
It shows that our method, especially without artificial friction, can
approximate the look of reference FLIP simulation qualitatively without the
noise on FLIP surface reconstructed from particles. The reference FLIP
simulation at half of the resolution of the hybrid FLIP already costs 120 Gb
memory with 1.7 billion particles. Memory consumption prevents us to compute
the reference at the same resolution. A photograph of real boat wakes is in
Figure 12. In Figure 11e, we show the result of a height field approach with
correct dispersion relationship (Canabal et al., 2016).

0.20.61.01.41.8610141822Froude NumberWake AngleOur MethodReal Data Figure 13.
Wake angle as a function of the hull Froude number $Fr=U/\sqrt{gL}$, where $U$
is the flow velocity, $g$ is the gravity, and $L$ is the length of the boat.
Our results match the measurements based on airborne images from (Rabaud and
Moisy, 2013).

Next we analyze the wake angles. The wakes behind a ship/boat are usually
referred to as the Kelvin wakes, which are characterized by a fixed angle
around $20^{\circ}$ and feather like pattern. The fixed angle is due to the
dispersion relationship of deep water waves of all wavelengths. However in
reality, Rabaud and Moisy (2013) reported the transition of fixed Kelvin angle
to increasingly smaller Mach angles when the ship/boat moves faster. They
believe ”the finite size of the disturbance (ship) explains this transition”.
Our hybrid approach reproduces their discoveries, matching their theory as
well as their real world measurements from airborne images. In Figure 13 we
can observe the transition with increasing speed.

Figure 14. FLIP-BEM transition boundary in the hybrid simulation in Figure
10d. The transition zone is visible in the first zoomed in figure, since
particles above the BEM mesh are removed in the coupling algorithm to match
the BEM free surface. Particles and triangles are colored by velocity.

In the transition zone between FLIP and BEM we fill and remove FLIP particles
to track the BEM free surface. This is visualized in Figure 14.

\begin{overpic}[width=203.80193pt]{graphics_final/overturning_hybrid_ref392.png}
\put(1.0,54.0){(a)} \put(51.0,54.0){(b)} \put(49.0,50.5){5 m/s}
\put(49.0,39.0){0 m/s} \end{overpic} Figure 15.  Extending overturning waves
behind a wedge with (a) particles and the BEM mesh in our hybrid method and
(b) particles in the reference FLIP simulation. Only the top layer of
particles are exported to save IO. The domain size is $20\,\textrm{m}\times
10\,\textrm{m}\times{45}\,\textrm{m}$. $\Delta
x_{\text{BEM}}=0.14\,\textrm{m},\Delta x_{\text{FLIP}}=0.07\,\textrm{m}$

Our method is capable of extending overturning waves behind the obstacle.
However, the BEM is not the best choice for such turbulent waves because it
handles complex topology changes less efficiently compared to the FLIP method.
In addition, the BEM is not designed for capturing vortices near the surface.
Hence extending overturning waves are not the recommended usage of our
approach. Nevertheless, we show the overturning wave results in Figure 15.

Figure 16. Simulation of a speedboat following a serpentine movement. The BEM
allows for the natural extension of the waves. Particles are rendered inside
the FLIP domain. The BEM mesh subtracting the FLIP domain patch are rendered
on the outside.

In a more complex example, the boat moves like a sine wave. The FLIP domain
size is $32\,\textrm{m}\times 3.5\,\textrm{m}\times 85\,\textrm{m}$. The
results are shown in Figure 16.

###  5.6. Multiple Moving Domains and Rendering

\begin{overpic}[width=433.62pt]{graphics_final/composed_triple_boat_passes.png}
\put(0.5,59.0){(a)} \put(50.5,59.0){(b)} \put(0.5,38.7){(c)}
\put(50.5,38.7){(d)} \put(0.5,18.4){(e)} \put(50.5,18.4){(f)} \end{overpic}
Figure 17. Simulation of speedboats riding in a circle computed with our
hybrid method. (a) Raw BEM mesh without subdivision colored by velocity. (b)
Rendered BEM mesh and exterior ocean with FFT waves. (c) Raw top layer FLIP
particles colored by velocity. (d) Rendered FLIP mesh with FFT waves. (e)
Composed image of b and d. (f) Composed image of a secondary whitewater
simulation. The water surface is perturbed by the previous FLIP domain and
affects the following FLIP domain.

In this example, we show that the method handles multiple moving FLIP domains
interacting with each other through BEM meshes. As shown in Figure 1, and 17,
three FLIP domains follow the cycling motion of three boats. The boat at the
front excites the water surface. When a FLIP domain in the back sweeps through
the perturbed region, newly generated particles follow the perturbed water
surface and its motion. There are three $50\,\textrm{m}\times
3.5\,\textrm{m}\times 50\,\textrm{m}$ FLIP domains in a cylinder tank with 100
m radius and 20 m depth. The track radius is 40 m while speedboats travel at
10 m/s. In principle, the FLIP domains can overlap, because they only interact
with each other through the BEM mesh. However, we avoid doing so because it
can introduce extra difficulty in rendering overlapping surfaces.

We use this scene as an example to demonstrate a simplified production
pipeline. Our hybrid method allows secondary whitewater simulation to enhance
the final look. This is possible because the BEM mesh provides surface
velocity field. We convert the BEM mesh to SDF with OpenVDB, and create narrow
band volumetric velocity field around the mesh according to the closest point
velocity on the BEM mesh for each voxel. The BEM SDF and velocity field are
replaced by FLIP SDFs and velocity fields inside FLIP domains. Such combined
velocity field and SDF are used for the off-the-shelf whitewater simulation
included in Houdini (Side Effects Software, 2021).

The final image is composed of different render passes. Initially we have a
BEM mesh (Figure17a). We smooth it and flatten the edge near BEM boundary to
match a flat ocean surface. The intersection of BEM mesh and a large cylinder
is rendered together with the exterior flat plane subtracting the same
cylinder to create seamless ocean. Next, FFT waves (Tessendorf, 2001) are
added at low velocity regions. This forms the first pass (Figure 17b).

We only output the surface layer of particles (Figure 17c) and the other
particles are exported as a SDF at simulation resolution to save disk
consumption. The surface particles are turned into a higher resolution SDF and
merged with the rest simulation resolution SDF. The combined SDF are turned
into polygons. We extract the FLIP surface by taking the intersection of three
FLIP polygon meshes and three boxes slightly smaller than the FLIP domains. In
order to reduce artifacts on the FLIP-BEM boundary, near the intersecting box,
we pull the intersected FLIP mesh vertices towards the nearest point on
smoothed BEM mesh. The closer they are to the intersecting box, the more they
cling to the smoothed BEM mesh. These flattened FLIP meshes are augmented with
FFT waves to render the high resolution liquid mesh (Figure 17d). When
rendering Figure 17d, the smoothed BEM mesh subtracting the three boxes
participates in light transport, but are not rendered (Houdini phantom).

Figure 17d is put on top of Figure 17b to get the rendered liquid surface in
Figure 17e. Finally the whitewater pass are added in Figure 17f.

##  6\. Runtime Performance

Table 2. Runtime comparison of our FLIP and Houdini (Side Effects Software,
2021), which implements narrow band FLIP (NB) (Ferstl et al., 2016) and
adaptive octree solver for pressure (Octree) (Goldade et al., 2019). The
comparison is on the same machine A in Table 1 with three substeps/frame. The
first two rows corresponds to the FLIP region in Figure 10h. The last two rows
corresponds to brute force simulation of high resolution ($\Delta
x_{\text{FLIP}}=0.07$ m) shallow ocean in Figure 10g ($\Delta
x_{\text{FLIP}}=0.14$ m). Our FLIP outperforms Houdini even with stricter
accuracy ($\varepsilon$) for the pressure solver.

Method | NB | Octree | DOF | Liquid Voxels | $\varepsilon$ | t/frame
---|---|---|---|---|---|---
Our FLIP | NO | NO | 12.2M | $285\times 50\times 857$ | $10^{-7}$ | 28s
Houdini | YES | YES | 2.4M | $285\times 50\times 857$ | $10^{-4}$ | 43s
YES | NO | 12.2M | $285\times 50\times 857$ | $10^{-4}$ | 90s
NO | NO | 12.2M | $285\times 50\times 857$ | $10^{-4}$ | 105s
Our FLIP | NO | NO | 293M | $1714\times 50\times 3429$ | $10^{-7}$ | 576s
Houdini | YES | YES | 50M | $1714\times 50\times 3429$ | $10^{-4}$ | 1077s

One of the motivations of our work is the limited runtime performance of
existing methods when simulating a large body of water. In this section we
compare to several other possible methods in terms of runtime performance. The
runtime performances of some methods are only available for different scenes
and machine configurations. Therefore we scale these runtimes linearly for
rough estimation.

First we compare variations of FLIP solvers. We compare our FLIP with narrow
band FLIP (Ferstl et al., 2016) and adaptive octree pressure solver (Goldade
et al., 2019) provided by Houdini (Side Effects Software, 2021) directly on
the same machine. The comparison is summarized in Table 2. First we compare
the FLIP part in our hybrid simulation of wake pattern (Figure 10). Our FLIP
implementation is faster than the two improvements combined. Furthermore we
try to simulate the whole ocean with brute force. However, neither Houdini nor
our FLIP can fit into the memory. Therefore we simulate the first three sub-
steps of the shallow ocean in Figure 10 but at high resolution, which amounts
to extending the FLIP region horizontally to fit the BEM region. Tall-cells
methods (Irving et al., 2006; Chentanez and Müller, 2011) are slower than
uniform methods for the same degrees of freedom (DOF). Even if they run as
fast as the uniform methods with the same DOF, given the runtime in Table 2, a
5-cell-thick simulation of the large domain with tall-cells methods would take
58s , still slower than our hybid method (42s in Table 1).

Next we roughly estimate the performance of a recent octree simulator (Ando
and Batty, 2020). It provides not only adaptivity inside the liquid volume,
but also on the liquid surface, making it very suitable for both deep and wide
liquid simulation. Our hybrid framework is comparable to its high performance.
They reported 20s/substep for a water tank scene with equivalently 33.3M
uniform liquid voxels (or 266M particles) on a 8-core Intel CPU at 3.6 GHz.

We estimate the run time if we use such octree simulator for the FLIP part of
our straight boat wake scene in Figure 10. Taking into account of CPU
frequency and core numbers, their simulator would take
$20\,\textrm{s/substep}\times 8\,\textrm{cores}/24\,\textrm{cores}\times
3.6\,\textrm{GHz}/3.0\,\textrm{GHz}\times
99\textrm{M}/266\textrm{M}=3.0\,\textrm{s/substep}$. Hence compared to our 9.3
s/substep FLIP performance in Table 1, it is three times faster. The BEM
domain is 24 times larger and 7 times deeper. It is possible to simulate this
using the octree method because these additional volumes would be filled with
coarse grids. However, whether it can preserve wakes and simulate the whole
domain efficiently cannot be concluded with our estimation. In their airplane
scene, the wakes seem to dissipate quickly.

##  7\. Discussion and Future Work

In this paper, we hybridize a local high-resolution FLIP solver and a surface-
based BEM solver for multiscale ocean modeling. This method allows us to
naturally extend the FLIP domain to a larger extent as shown in the crown
splash experiment (Figure 9). The accuracy is demonstrated by comparing the
boat wake pattern to the brute-force FLIP results and a real example (Figure
10, 11 13). Through the BEM mesh, multiple volumetric domains can interact
with each other and collectively generate complex patterns in a large body of
water (Figure 17). Thanks to the surface representation of liquid and adaptive
triangle sizes of the BEM, the BEM mesh can cheaply expand both horizontally
and vertically.

Our method is suitable for open water simulation with both water splashes near
the main object, and waves that travel long distances. The most time consuming
part of the simulation is the volumetric solver near the main object, even
with our optimized implementation. The extra time induced by the BEM
simulation in the calm region and coupling is usually no more than $30\%$. The
FLIP part consumes about $90\%$ of the memory, up to 120 Gb for the 1.7 billon
particles scene. It is up to the user to decide the resolution and size of the
FLIP region in the BEM water body depending on the machine capability. The BEM
is limited by the acceleration grid for collision detection. Currently, it
uses uniform grids with a cell size equal to the smallest triangle edge
length. Hence, the bounding box of the BEM cannot be too large.

The performance of the BEM is the main drawback. It is limited by the serial
surface tracking program LosTopos (Da et al., 2014), which takes quite a long
time to resolve collisions and topology changes, etc. The BEM involves dense
matrix-vector multiplication as a result of integral equation systems. In
principle, the fast multipole method (FMM) (Greengard and Rokhlin, 1987) can
solve the matrix-vector multiplication problem, but for the problem sizes in
our experiments, the naive summation on the GPU is faster given the overhead
of the FMM. In our most challenging triple boat scene, the BEM simulation
takes up to 2 minutes per frame for 170K triangles. About half of the time is
spent on remeshing. Therefore, the smart usage of the BEM is for calm waves,
rather than splashes and overturning waves. When the main object is placed in
a stormy ocean, the propagated waves are negligible. The prescribed waves with
one-way interaction are more efficient. The BEM does not interacts with
beaches, cliffs either, because of the irregularity of solids unsupported by
LosTopos. Another issue is the wake amplitudes in Figure 10. We alleviate this
issue by sampling below the interface. We believe this problem is due to the
Helmholtz decomposition which interprets slow velocity as resistance. The
resistance causes the elevation of the BEM mesh ahead of the boat. When the
detailed water waves in the FLIP domain enter the BEM domain, most of the
details will be lost because the BEM mesh has lower resolution, dedicated for
large wave features.

The coupling workflow may look heuristic compared to a strong coupling
framework that bridges a dense integral equation system in BEM and a sparse
equation system in FLIP. However, at this stage we are satisfied with its
results, considering potentially larger implementation effort of a strong
coupling system, and corresponding lower efficiency of a mixed integral and
differential equation system.

For the problem of extending waves with correct dispersion relationships, we
believe the dispersion kernel method (Canabal et al., 2016) shows great
potential. We look forward to seeing coupling of such height field methods
with correct dispersion relationships with 3D simulations. Canabal et al.
directly modified the height field, so the height field does not have
horizontal displacement or velocity information. The two-way coupling for
domain extension may not be trivial due to their inherently different physics.

It would be interesting to see our method applied to other types of hybrid
Lagrangian-Eulerian simulations such as IQ-MPM (Fang et al., 2020).

## Acknowledgements

This work has been supported by KAUST (individual baseline funding of the
Computational Sciences Group within the Visual Computing Center), the National
Science Foundation (NSF CAREER IIS-1943199), and CCF-1813624 and ECCS-2023780.

We thank Da Fang and Christopher Batty for releasing their source code and
Miguelangelo Rosario for providing the ship model. The valuable comments of
the anonymous reviewers are gratefully acknowledged.

## References

  * (1)
  * Aanjaneya et al. (2017) Mridul Aanjaneya, Ming Gao, Haixiang Liu, Christopher Batty, and Eftychios Sifakis. 2017.  Power Diagrams and Sparse Paged Grids for High Resolution Adaptive Liquids.  _ACM Trans. Graph._ 36, 4, Article 140 (July 2017), 12 pages.
  * Ando and Batty (2020) Ryoichi Ando and Christopher Batty. 2020.  A Practical Octree Liquid Simulator with Adaptive Surface Resolution.  _ACM Trans. Graph._ 39, 4, Article 32 (July 2020), 17 pages.
  * Ando et al. (2013) Ryoichi Ando, Nils Thuerey, and Chris Wojtan. 2013\.  Highly Adaptive Liquid Simulations on Tetrahedral Meshes.  _ACM Trans. Graph._ 32, 4, Article 103 (July 2013), 10 pages.
  * Batty (2018) Christopher Batty. 2018\.  Fluid3D.  https://github.com/christopherbatty/Fluid3D.
  * Batty et al. (2007) Christopher Batty, Florence Bertails, and Robert Bridson. 2007\.  A Fast Variational Framework for Accurate Solid-Fluid Coupling.  _ACM Trans. Graph._ 26, 3 (July 2007), 100–es.
  * Bhattacharya et al. (2012) Haimasree Bhattacharya, Michael Bang Nielsen, and Robert Bridson. 2012.  Steady State Stokes Flow Interpolation for Fluid Control.. In _Eurographics (Short Papers)_. Citeseer, 57–60.
  * Bojsen-Hansen and Wojtan (2016) Morten Bojsen-Hansen and Chris Wojtan. 2016.  Generalized Non-Reflecting Boundaries for Fluid Re-Simulation.  _ACM Trans. Graph._ 35, 4, Article 96 (July 2016), 7 pages.
  * Brackbill and Ruppel (1986) Jeremiah U. Brackbill and Hans M. Ruppel. 1986.  FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions.  _J. Comput. Phys._ 65, 2 (1986), 314–343.
  * Bridson and Müller-Fischer (2007) Robert Bridson and Matthias Müller-Fischer. 2007\.  Fluid Simulation: SIGGRAPH 2007 Course Notes.  In _ACM SIGGRAPH 2007 Courses_. Association for Computing Machinery, New York, NY, USA, 1–81.
  * Canabal et al. (2016) José A. Canabal, David Miraut, Nils Thuerey, Theodore Kim, Javier Portilla, and Miguel A. Otaduy. 2016\.  Dispersion Kernels for Water Wave Simulation.  _ACM Trans. Graph._ 35, 6, Article 202 (Nov. 2016), 10 pages.
  * Chen et al. (2018) Xi Chen, Ren chuan Zhu, Wen jun Zhou, and Ji Zhao. 2018\.  A 3D multi-domain high order boundary element method to evaluate time domain motions and added resistance of ship in waves.  _Ocean Engineering_ 159 (2018), 112 – 128.
  * Chentanez et al. (2015) Nuttapong Chentanez, Matthias Müller, and Tae-Yong Kim. 2015\.  Coupling 3D eulerian, heightfield and particle methods for interactive simulation of large scale liquid phenomena.  _IEEE Transactions on Visualization and Computer Graph (TVCG)_ 21, 10 (2015), 1116–1128.
  * Chentanez and Müller (2011) Nuttapong Chentanez and Matthias Müller. 2011.  Real-Time Eulerian Water Simulation Using a Restricted Tall Cell Grid.  _ACM Trans. Graph._ 30, 4, Article 82 (July 2011), 10 pages.
  * Colagrossi and Landrini (2003) Andrea Colagrossi and Maurizio Landrini. 2003.  Numerical Simulation of Interfacial Flows by Smoothed Particle Hydrodynamics.  _J. Comput. Phys._ 191, 2 (Nov. 2003), 448–475.
  * Da et al. (2014) Fang Da, Christopher Batty, and Eitan Grinspun. 2014\.  Multimaterial Mesh-Based Surface Tracking.  _ACM Trans. Graph._ 33, 4, Article 112 (July 2014), 11 pages.
  * Da et al. (2016a) Fang Da, David Hahn, Christopher Batty, Chris Wojtan, and Eitan Grinspun. 2016a.  Surface-Only Liquids.  _ACM Trans. Graph._ 35, 4, Article 78 (July 2016), 12 pages.
  * Da et al. (2016b) Fang Da, David Hahn, Christopher Batty, Chris Wojtan, and Eitan Grinspun. 2016b.  Surface-Only Liquids: Source Code.  http://www.cs.columbia.edu/cg/surfaceliquids/code/.
  * Edmont (2009) Edmont. 2009.  _Fjordn surface wave boat_.  https://commons.wikimedia.org/wiki/File:Fjordn_surface_wave_boat.jpg
  * English et al. (2013) R. Elliot English, Linhai Qiu, Yue Yu, and Ronald Fedkiw. 2013\.  Chimera Grids for Water Simulation. In _Proceedings of the 12th ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Anaheim, California) _(SCA ’13)_. Association for Computing Machinery, New York, NY, USA, 85–94.
  * Fang et al. (2020) Yu Fang, Ziyin Qu, Minchen Li, Xinxin Zhang, Yixin Zhu, Mridul Aanjaneya, and Chenfanfu Jiang. 2020.  IQ-MPM: An Interface Quadrature Material Point Method for Non-Sticky Strongly Two-Way Coupled Nonlinear Solids and Fluids.  _ACM Trans. Graph._ 39, 4, Article 51 (July 2020), 16 pages.
  * Ferstl et al. (2016) Florian Ferstl, Ryoichi Ando, Chris Wojtan, Rüdiger Westermann, and Nils Thuerey. 2016\.  Narrow Band FLIP for Liquid Simulations. In _Proceedings of the 37th Annual Conference of the European Association for Computer Graphics_ (Lisbon, Portugal) _(EG ’16)_. Eurographics Association, Goslar, DEU, 225–232.
  * Fu et al. (2017) Chuyuan Fu, Qi Guo, Theodore Gast, Chenfanfu Jiang, and Joseph Teran. 2017.  A Polynomial Particle-in-Cell Method.  _ACM Trans. Graph._ 36, 6, Article 222 (Nov. 2017), 12 pages.
  * Gibou et al. (2002) Frederic Gibou, Ronald P Fedkiw, Li-Tien Cheng, and Myungjoo Kang. 2002.  A second-order-accurate symmetric discretization of the Poisson equation on irregular domains.  _J. Comput. Phys._ 176, 1 (2002), 205–227.
  * Golas et al. (2012) Abhinav Golas, Rahul Narain, Jason Sewall, Pavel Krajcevski, Pradeep Dubey, and Ming Lin. 2012\.  Large-Scale Fluid Simulation Using Velocity-Vorticity Domain Decomposition.  _ACM Trans. Graph._ 31, 6, Article 148 (Nov. 2012), 9 pages.
  * Goldade et al. (2019) Ryan Goldade, Yipeng Wang, Mridul Aanjaneya, and Christopher Batty. 2019.  An Adaptive Variational Finite Difference Framework for Efficient Symmetric Octree Viscosity.  _ACM Trans. Graph._ 38, 4, Article Article 94 (July 2019), 14 pages.
  * Greengard and Rokhlin (1987) Leslie Greengard and Vladimir Rokhlin. 1987.  A fast algorithm for particle simulations.  _Journal of computational physics_ 73, 2 (1987), 325–348.
  * Grilli et al. (1989) Stephan T. Grilli, J. Skourup, and I.A. Svendsen. 1989\.  An efficient boundary element method for nonlinear water waves.  _Engineering Analysis with Boundary Elements_ 6, 2 (1989), 97 – 107.
  * Hahn and Wojtan (2016) David Hahn and Chris Wojtan. 2016.  Fast Approximations for Boundary Element Based Brittle Fracture Simulation.  _ACM Trans. Graph._ 35, 4, Article 104 (July 2016), 11 pages.
  * Harlow (1964) Francis H Harlow. 1964\.  The particle-in-cell computing method for fluid dynamics.  _Methods Comput. Phys._ 3 (1964), 319–343.
  * Harlow and Welch (1965) Francis H. Harlow and J. Eddie Welch. 1965.  Numerical Calculation of Time‐Dependent Viscous Incompressible Flow of Fluid with Free Surface.  _Physics of Fluids_ 8 (1965), 2182–2189.
  * Hong (2009) Woo-Suck Hong. 2009\.  _An Adaptive Sampling Approach to Incompressible Particle-Based Fluid_.  Texas A & M University, USA.  AAI3370710.
  * Huang et al. (2019) Libo Huang, Torsten Hädrich, and Dominik L. Michels. 2019\.  On the Accurate Large-scale Simulation of Ferrofluids.  _ACM Trans. Graph._ 38, 4, Article 93 (July 2019), 15 pages.
  * Huang and Michels (2020) Libo Huang and Dominik L. Michels. 2020.  Surface-Only Ferrofluids.  _ACM Trans. Graph._ 39, 6, Article 174 (Nov. 2020), 17 pages.
  * Ibayashi et al. (2020) H. Ibayashi, C. Wojtan, N. Thuerey, T. Igarashi, and R. Ando. 2020.  Simulating Liquids on Dynamically Warping Grids.  _IEEE Transactions on Visualization and Computer Graphics_ 26, 06 (jun 2020), 2288–2302.
  * Ihmsen et al. (2012) Markus Ihmsen, Nadir Akinci, Gizem Akinci, and Matthias Teschner. 2012.  Unified spray, foam and air bubbles for particle-based fluids.  _The Visual Computer_ 28, 6 (2012), 669–677.
  * Ihmsen et al. (2014) Markus Ihmsen, Jens Cornelis, Barbara Solenthaler, Christopher Horvath, and Matthias Teschner. 2014.  Implicit Incompressible SPH.  _IEEE Transactions on Visualization and Computer Graphics_ 20, 3 (2014), 426–435.
  * Irving et al. (2006) Geoffrey Irving, Eran Guendelman, Frank Losasso, and Ronald Fedkiw. 2006.  Efficient Simulation of Large Bodies of Water by Coupling Two and Three Dimensional Techniques.  _ACM Trans. Graph._ 25, 3, 805–811.
  * James and Pai (1999) Doug L. James and Dinesh K. Pai. 1999.  ArtDefo: Accurate Real Time Deformable Objects. In _Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques_ _(SIGGRAPH ’99)_. ACM Press/Addison-Wesley Publishing Co., USA, 65–72.
  * Jeschke et al. (2020) Stefan Jeschke, Christian Hafner, Nuttapong Chentanez, Miles Macklin, Matthias Müller-Fischer, and Christopher Wojtan. 2020.  Making Procedural Water Waves Boundary-aware.  _Computer Graphics Forum_ 39, 8 (2020), 47–54.
  * Jeschke et al. (2018) Stefan Jeschke, Tomáš Skřivan, Matthias Müller-Fischer, Nuttapong Chentanez, Miles Macklin, and Chris Wojtan. 2018\.  Water Surface Wavelets.  _ACM Trans. Graph._ 37, 4, Article 94 (July 2018), 13 pages.
  * Jeschke and Wojtan (2015) Stefan Jeschke and Chris Wojtan. 2015.  Water Wave Animation via Wavefront Parameter Interpolation.  _ACM Trans. Graph._ 34, 3, Article 27 (May 2015), 14 pages.
  * Jeschke and Wojtan (2017) Stefan Jeschke and Chris Wojtan. 2017.  Water Wave Packets.  _ACM Trans. Graph._ 36, 4, Article 103 (July 2017), 12 pages.
  * Jiang et al. (2015) Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin. 2015.  The Affine Particle-in-Cell Method.  _ACM Trans. Graph._ 34, 4, Article 51 (July 2015), 10 pages.
  * Keeler and Bridson (2015) Todd Keeler and Robert Bridson. 2015.  Ocean Waves Animation Using Boundary Integral Equations and Explicit Mesh Tracking. In _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Copenhagen, Denmark) _(SCA ’14)_. Eurographics Association, Goslar, DEU, 11–19.
  * Kinnas and Hsin (1992) Spyros A. Kinnas and Ching-Yeh Hsin. 1992.  Boundary element method for the analysis of the unsteady flow aroundextreme propeller geometries.  _AIAA Journal_ 30, 3 (1992), 688–696.
  * Longuet-Higgins and Cokelet (1976) Michael Selwyn Longuet-Higgins and ED Cokelet. 1976.  The deformation of steep surface waves on water-I. A numerical method of computation.  _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ 350, 1660 (1976), 1–26.
  * Losasso et al. (2004) Frank Losasso, Frédéric Gibou, and Ronald Fedkiw. 2004.  Simulating Water and Smoke with an Octree Data Structure.  _ACM Trans. Graph._ 23, 3 (Aug. 2004), 457–462.
  * Losasso et al. (2008) Frank Losasso, Jerry Talton, Nipun Kwatra, and Ronald Fedkiw. 2008.  Two-Way Coupled SPH and Particle Level Set Fluid Simulation.  _IEEE Transactions on Visualization and Computer Graphics_ 14, 4 (2008), 797–804.
  * Museth (2013) Ken Museth. 2013\.  VDB: High-Resolution Sparse Volumes with Dynamic Topology.  _ACM Trans. Graph._ 32, 3, Article 27 (July 2013), 22 pages.
  * Ng et al. (2009) Yen Ting Ng, Chohong Min, and Frédéric Gibou. 2009\.  An Efficient Fluid-Solid Coupling Algorithm for Single-Phase Flows.  _J. Comput. Phys._ 228, 23 (Dec. 2009), 8807–8829.
  * Nielsen and Bridson (2011) Michael B. Nielsen and Robert Bridson. 2011.  Guide Shapes for High Resolution Naturalistic Liquid Simulation.  _ACM Trans. Graph._ 30, 4, Article 83, 8 pages.
  * Nielsen and Bridson (2016) Michael B. Nielsen and Robert Bridson. 2016.  Spatially Adaptive FLIP Fluid Simulations in Bifrost. In _ACM SIGGRAPH 2016 Talks_ (Anaheim, California) _(SIGGRAPH ’16)_. Association for Computing Machinery, New York, NY, USA, Article 41, 2 pages.
  * Nielsen et al. (2013) Michael B. Nielsen, Andreas Söderström, and Robert Bridson. 2013.  Synthesizing Waves from Animated Height Fields.  _ACM Trans. Graph._ 32, 1, Article 2 (Feb. 2013), 9 pages.
  * Nielsen et al. (2017) Michael B. Nielsen, Konstantinos Stamatelos, Adrian Graham, Marcus Nordenstam, and Robert Bridson. 2017\.  Localized Guided Liquid Simulations in Bifrost. In _ACM SIGGRAPH 2017 Talks_ (Los Angeles, California) _(SIGGRAPH ’17)_. Association for Computing Machinery, New York, NY, USA, Article 44, 2 pages.
  * Rabaud and Moisy (2013) Marc Rabaud and Frédéric Moisy. 2013.  Ship wakes: Kelvin or Mach angle?  _Physical review letters_ 110, 21 (2013), 214503\.
  * Schreck et al. (2019) Camille Schreck, Christian Hafner, and Chris Wojtan. 2019\.  Fundamental Solutions for Water Wave Animation.  _ACM Trans. Graph._ 38, 4, Article 130 (July 2019), 14 pages.
  * Setaluri et al. (2014) Rajsekhar Setaluri, Mridul Aanjaneya, Sean Bauer, and Eftychios Sifakis. 2014.  SPGrid: A Sparse Paged Grid Structure Applied to Adaptive Smoke Simulation.  _ACM Trans. Graph._ 33, 6, Article 205 (Nov. 2014), 12 pages.
  * Side Effects Software (2021) Side Effects Software. 2021\.  Houdini.  (2021).
  * Söderström et al. (2010) Andreas Söderström, Matts Karlsson, and Ken Museth. 2010\.  A PML-Based Nonreflective Boundary for Free Surface Fluid Animation.  _ACM Trans. Graph._ 29, 5, Article 136 (Nov. 2010), 17 pages.
  * Solenthaler et al. (2011) Barbara Solenthaler, Peter Bucher, Nuttapong Chentanez, Matthias Müller, and Markus Gross. 2011\.  SPH Based Shallow Water Simulation. In _Workshop in Virtual Reality Interactions and Physical Simulation ”VRIPHYS” (2011)_. The Eurographics Association.
  * Solomon et al. (2017) Justin Solomon, Amir Vaxman, and David Bommes. 2017\.  Boundary Element Octahedral Fields in Volumes.  _ACM Trans. Graph._ 36, 3, Article 28 (May 2017), 16 pages.
  * Stüben (2001) Klaus Stüben. 2001\.  A review of algebraic multigrid.  _J. Comput. Appl. Math._ 128, 1 (2001), 281–309.  Numerical Analysis 2000. Vol. VII: Partial Differential Equations.
  * Stomakhin and Selle (2017) Alexey Stomakhin and Andrew Selle. 2017.  Fluxed Animated Boundary Method.  _ACM Trans. Graph._ 36, 4, Article 68 (July 2017), 8 pages.
  * Tessendorf (2001) Jerry Tessendorf. 2001\.  Simulating Ocean Water.  _SIGGRAPH’99 Course Note_ (01 2001).
  * Thuerey et al. (2007) Nils Thuerey, Matthias Müller-Fischer, Simon Schirm, and Markus Gross. 2007.  Real-time Breaking Waves for Shallow Water Simulations. In _15th Pacific Conference on Computer Graphics and Applications (PG’07)_. 39–46.
  * Thuerey et al. (2006) Nils Thuerey, Ulrich Rüde, and Marc Stamminger. 2006\.  Animation of Open Water Phenomena with coupled Shallow Water and Free Surface Simulations. In _ACM SIGGRAPH / Eurographics Symposium on Computer Animation_ , Marie-Paule Cani and James O’Brien (Eds.). The Eurographics Association.
  * Weber et al. (2015) Daniel Weber, Johannes Mueller-Roemer, André Stork, and Dieter Fellner. 2015.  A Cut-Cell Geometric Multigrid Poisson Solver for Fluid Simulation.  _Computer Graphics Forum_ 34, 2 (2015), 481–491.
  * Winchenbach et al. (2017) Rene Winchenbach, Hendrik Hochstetter, and Andreas Kolb. 2017\.  Infinite Continuous Adaptivity for Incompressible SPH.  _ACM Trans. Graph._ 36, 4, Article 102 (July 2017), 10 pages.
  * Wu and Taylor (2003) Guoxiong Wu and R. Eatock Taylor. 2003.  The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies.  _Ocean Engineering_ 30, 3 (2003), 387 – 400.
  * Yue et al. (2018) Yonghao Yue, Breannan Smith, Peter Yichen Chen, Maytee Chantharayukhonthorn, Ken Kamrin, and Eitan Grinspun. 2018.  Hybrid Grains: Adaptive Coupling of Discrete and Continuum Simulations of Granular Media.  _ACM Trans. Graph._ 37, 6, Article 283, 19 pages.
  * Yuksel et al. (2007) Cem Yuksel, Donald H. House, and John Keyser. 2007\.  Wave Particles. In _ACM SIGGRAPH 2007 Papers_ (San Diego, California) _(SIGGRAPH ’07)_. Association for Computing Machinery, New York, NY, USA, 99–es.
  * Zhang and Liu (2015) Aman Zhang and Yun-Long Liu. 2015.  Improved three-dimensional bubble dynamics model based on boundary element method.  _J. Comput. Phys._ 294 (2015), 208 – 223.
  * Zheng and James (2010) Changxi Zheng and Doug L. James. 2010.  Rigid-Body Fracture Sound with Precomputed Soundbanks.  In _ACM SIGGRAPH 2010 Papers_. Association for Computing Machinery, New York, NY, USA.
  * Zhu et al. (2013) Bo Zhu, Wenlong Lu, Matthew Cong, Byungmoon Kim, and Ronald Fedkiw. 2013.  A New Grid Structure for Domain Extension.  _ACM Trans. Graph._ 32, 4, Article 63 (July 2013), 12 pages.
  * Zhu and Bridson (2005) Yongning Zhu and Robert Bridson. 2005.  Animating Sand as a Fluid.  _ACM Trans. Graph._ 24, 3 (July 2005), 965–972.