Water surface wavelets

Works in both physics and computer graphics employ numerous theoretical assumptions to the Navier-Stokes equations in order to reduce the complexity and make them analytically tractable. Some common assumptions in computer graphics are deep water, potential flow, small amplitudes, and periodic boundary conditions. These approaches take advantage of analytical solutions to the Navier-Stokes equations on a 2D height field in order to express the motion of the ocean in the form of sines and cosines, while sacrificing the ability to simulate arbitrary fluid motion [Hinsinger et al. 2002; Horvath 2015; Mastin et al. 1987; Tessendorf 2004b]. Jeschke & Wojtan [2015] simulate the propagation of wavefronts throughout a static environment, in order to extend spectrum-based approaches to handle complex boundaries. However, the method is limited to pre-computation, and does not allow interactions with moving obstacles.

We call this technique of discretizing a Fourier-transform the “spectrum-based” approach. These methods exhibit theoretically unlimited visual detail, in the sense that they can animate arbitrarily high frequency waves without impacting the method’s accuracy or stability. Similarly, fast wave speeds are trivial to simulate, because the motion is independent of the time step size. However, because the derivation of these methods makes several assumptions about the underlying flow, they tend to have limitations like the inability to realistically interact with complex boundaries.

A great way to get around the limitations of spectrum-based methods is to directly simulate a two-dimensional version of the Navier-Stokes equations using numerical algorithms. Some approaches discretize simplified wave equations [Kass and Miller 1990; Thuerery et al. 2010; Yang et al. 2016] or thin plate equations [Yu et al. 2012]; these simplified equations are easier to implement, but the resulting behavior is fundamentally different from actual water waves. The work of Tessendorf [2004a; 2014], which discretizes a linearized Bernoulli equation, exhibits more realistic wave dispersion but still falls short of physically correct motion (as discussed in [Canabal et al. 2016]). Some researchers have also introduced direct numerical simulators based on the lattice Boltzmann method (LBM) [Geist et al. 2010], which requires a careful tuning of the LBM collision matrix to yield realistic wave speeds. Convolution-based approaches [Loviscach 2002; Ottosson 2011] aim to achieve the correct dispersion relation, but they must cope with the practical difficulties of large kernels occupying the entire simulation domain. Recently, Canabal et al. [2016] overcame these difficulties with a combination of pyramid filters and shadowed convolution operations.

Each of these methods tend to be much more flexible than spectrum-based methods, because they make minimal assumptions about the environment. For example, spectrum-based methods have difficulties with simulating obstacle interactions because their derivation assumes periodic boundaries, while direct simulation approaches have no such limitation. On the other hand, these methods must simulate wave propagation by iterating local kernel operations instead of simply plugging a time parameter into a cosine function. All of these methods discretize wave heights and momentum directly on an Eulerian grid or mesh, and, as a consequence, the resolution of the grid is directly tied to the amount of visible detail in the waves. Nyquist’s theorem requires that the grid must be fine enough to resolve the highest frequency in the heightfield; otherwise aliasing and instability may occur. Similarly, the stability of an explicitly integrated wave simulation is also intimately related to the grid spacing by the CFL condition. This stability problem can be overcome by implicit integration, at the expense of additional computation and a more complex and less GPU-friendly implementation.

Fully three-dimensional techniques for liquid simulation are outside the scope of this work. We recommend that interested readers consult the text by Bridson [2015].

Our method combines the flexibility of numerical approaches with the stability and visual detail of spectrum-based approaches, but we are not the first to do so. Yuksel et al. [2007] proposes wave particles which represent a local wave crest and move with a pre-determined wave speed $c$. We can view this approach as a local spectrum-based method if we use many different particle sizes and set $c$ equal to the analytical phase speed. Taking this further, Jeschke & Wojtan [2017] introduced wave packets which propagate at the theoretical group speed and contain a small train of waves traveling at the theoretical phase speed. These approaches inherit some advantages of spectrum-based methods, like the numerical stability and theoretically accurate wave speeds. At the same time, they avoid the complications with spectrum-based waves by breaking the global cosine waves into a train of shorter wave components that are free to interact with obstacles. Our method provides similar advantages, but it is Eulerian rather than Lagrangian (its degrees of freedom are associated with regions of space rather than the wave motion itself). Consequently, our method maps to GPU hardware more easily, the computational complexity is constant as it does

not vary with a number of particles, and it trivially interfaces with texture maps for easy artistic control.

Previous works have investigated liquid control by directly editing keyframes [McNamara et al. 2004; Shi and Yu 2005], by sculpting fluid interactively [Manteaux et al. 2016; Pan et al. 2013], by introducing guide forces on particles [Thurey et al. 2009] or meshes [Raveendran et al. 2012], by composing together convenient flow primitives [Chenney 2004], or by interpolating between simulations and existing motions [Raveendran et al. 2014; Thurey 2016]. Some works enhance existing animations with additional waves on surfaces [Angst et al. 2008; Kim et al. 2013] or near boundaries [Jeschke and Wojtan 2015, 2017]. Horvath et al. [Horvath 2015] explores the control of ocean spectra in great detail, and Nielsen et al. [Nielsen et al. 2013] investigates how to direct a spectrum-based method based on a scripted height field. We are unaware of any prior methods that locally direct the amplitudes of directional water waves, as we propose.

Current numerical methods for water wave simulation rely either on discretizations of partial differential equations (PDEs), or on spectrum-based methods. Methods based on discretized PDEs approximate some differential equation that describes how the wave height $\backslash eta(\backslash mathbf\{x\},\; t)$ evolves over time:

where the right hand side is defined by the particular wave model (shallow water equations, Bernoulli equation, etc.), and environmental interactions are encoded in the PDE’s boundary conditions. These discretizations sample $\backslash eta(\backslash mathbf\{x\},\; t)$ over space with a grid spacing equal to $\backslash Delta\; x$. In order to avoid aliasing and faithfully reproduce high-frequency details, The Nyquist-Shannon sampling theorem [Shannon 1949] requires that $\backslash Delta\; x$ is less than half of the shortest wavelength in $\backslash eta(\backslash mathbf\{x\},\; t)$. If $\backslash eta(\backslash mathbf\{x\},\; t)$ contains interesting high-frequency details, then the samples of $\backslash eta(\backslash mathbf\{x\},\; t)$ must be very close together, and thus $\backslash Delta\; x$ must be very small. Practically, decreasing $\backslash Delta\; x$ requires either more samples or a smaller simulation domain, so the Nyquist theorem effectively limits the visual detail that a simulator can produce. Highly detailed visuals can be obtained by decreasing $\backslash Delta\; x$ at the expense of vastly increased computation time and memory.

Spectrum-based methods for animating water waves [Tessendorf 2004b] remove these problems by avoiding a spatial discretization altogether. They rely on linear wave theory [Johnson 1997], which describes the wave height dynamics in terms of frequencies instead of partial derivatives:

Here, $\backslash eta\_c(\backslash mathbf\{x\},\; t)$ is a complex function that varies over two-dimensional space and time, and we can get the wave height by taking its real part, $\backslash eta(\backslash mathbf\{x\},\; t)\; =\; \backslash text\{Re\; \}\; \backslash eta\_c(\backslash mathbf\{x\},\; t)$. The wavevector $\backslash mathbf\{k\}$ is a two-dimensional frequency function, the wavenumber $k\; =\; |\backslash mathbf\{k\}|$ represents a scalar frequency, and $\backslash mathbf\{k\}\; =\; \backslash mathbf\{k\}/k$ is the wave direction. The exponential term in this equation represents a traveling wave, and $A(\backslash mathbf\{k\})$ represents its amplitude. The angular frequency

encodes the speed of each wave based on its wavenumber $k$, gravity $g$, and surface tension $\backslash sigma$. Spectrum-based methods compute the wave height by discretizing the integral over all of these two-dimensional waves, instead of discretizing a differential equation. This solution works perfectly for ideal situations with periodic domains, and no boundaries or interacting obstacles. However, these spectrum-based methods become impractical or impossible in less constrained scenarios.

To summarize, PDE-based methods for animating water waves excel in the simulation of low-frequency wavefunctions with complicated environmental interactions, while spectrum-based methods are ideal for simulating highly detailed wavefunctions undergoing simple motion without any boundary interactions. In the following section, we derive a hybrid discretization that relies on discretized PDEs to simulate low frequency motion and uses spectrum-based techniques for simulating the high frequencies. This strategy allows us to simulate wavelengths far shorter than $\backslash Delta\; x$ interacting with complex environments, while eliminating the previously mentioned problems associated with simulating highly detailed waves.

Our derivation relies on the Gabor wavelet transform, which effectively transforms Equation 2 (with amplitudes $A$ depending only on $\backslash mathbf\{k\}$ and independent of $\backslash mathbf\{x\}$ and $t$) to a similar one with amplitudes $\backslash mathcal\{A\}$ depending on $\backslash mathbf\{k\}$, $\backslash mathbf\{x\}$, and $t$:

Not only does this new spatially-varying amplitude $\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)$ allow more local control over the simulation compared to $A(\backslash mathbf\{k\})$, but we prove that $\backslash mathcal\{A\}$ is guaranteed to have lower frequency content than the wave height function $\backslash eta(\backslash mathbf\{x\},\; t)$. We will take advantage of this fact to forestall problems related to the Nyquist limit and create highly efficient and detailed wave simulations.

The Gabor transform [Gabor 1946] of the water height $\backslash eta\_c$ is

where the first exponential term in the integrand is a Gaussian centered at position $\backslash mathbf\{x\}$ with standard deviation $s$, and the second exponential term is a static wave with wavevector $\backslash mathbf\{k\}$. We can think of Equation 5 as an inner product between $\backslash eta\_c$ and a Gaussian wavelet,

so the number $\backslash zeta(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)$ tells how much the water height $\backslash eta\_c$ behaves like a wave with wave-vector $\backslash mathbf\{k\}$ in the vicinity of point $\backslash mathbf\{x\}$.

We can also invert the Gabor transform:

which reminds us of a Fourier transform with $\backslash zeta$ acting like an amplitude that now depends on $\backslash mathbf\{x\}$ and $t$ as well as $\backslash mathbf\{k\}$. We introduce the change of variables $\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)\; =\; \backslash zeta(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)\; e^\{i\backslash omega(\backslash mathbf\{k\})t\}$ to obtain an analogue to the dynamic wave evolution in Equation 2:

now with $\backslash mathcal\{A\}$ playing the role of an amplitude that varies over space and time as well as wavevector. By plugging Equation 2 into Equation 5, we obtain an equation for the time evolution of $\backslash mathcal\{A\}$:

Although this equation completely prescribes the evolution of $\backslash mathcal\{A\}$, its integral form makes it difficult to impose boundary conditions and generally inconvenient to discretize. We prefer to work with a differential equation instead, which we derive by taking the time derivative of Equation 8 and expressing it in terms of its spatial derivative. Linearizing the dispersion relation $\backslash omega(\backslash mathbf\{l\})\; \backslash approx\; \backslash omega(\backslash mathbf\{k\})\; +\; \backslash omega\text{'}(\backslash mathbf\{k\})\backslash mathbf\{k\}\; \backslash cdot\; (\backslash mathbf\{l\}\; -\; \backslash mathbf\{k\})$ gives us a first order partial differential equation for the evolution of $\backslash mathcal\{A\}$:

which tells us that our amplitude function $\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)$ gets advected in space in the direction $\backslash hat\{\backslash mathbf\{k\}\}$ with speed $\backslash omega\text{'}(\backslash mathbf\{k\})$. Note that this speed corresponds exactly to the group speed which transports water wave energy [Johnson 1997] and wave packets [Jeschke and Wojtan 2017]. We provide a more detailed derivation of this equation in Appendix A.

Equation 10 is subject to boundary conditions:

$\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)\; =\; \backslash mathcal\{A\}\_\{\backslash text\{ambient\}\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)$ on transmitting boundary

$\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)\; =\; \backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\}\_\{\backslash text\{reflect\}\},\; t)$ on reflecting boundary (11)

where $\backslash mathcal\{A\}\_\{\backslash text\{ambient\}\}$ is an amplitude function defined outside of the domain, and $\backslash mathbf\{k\}\_\{\backslash text\{reflect\}\}\; =\; \backslash mathbf\{k\}\; -\; 2(\backslash mathbf\{n\}\; \backslash cdot\; \backslash mathbf\{k\})\backslash mathbf\{n\}$ is the wavevector reflected off the boundary with a normal $\backslash mathbf\{n\}$. The reflecting boundary condition is based on how a planar wave reflects of a straight boundary.

Equations 7 and 10 are the key ingredients we need to simulate water waves. The actual interactive wave simulator advects $\backslash mathcal\{A\}$ through space each time step using Equation 10, and the renderer then uses Equation 7 to reconstruct the water surface where needed. In this approach, $\backslash mathcal\{A\}$ acts as the primary simulation variable, while $\backslash eta$ is used afterwards for reconstruction and visualization. ($\backslash mathcal\{A\}$ is updated by the simulation, and then we compute $\backslash eta$ where needed by summing together many different waves using $\backslash mathcal\{A\}$ as their amplitude.)

We note that the Gabor transform in Equation 5 is only used for the derivation, and there is no need to compute it. Similarly, the size of the Gaussian $s$ does not appear in Equations 7 and 10, so it is merely a theoretical construct that is useful for analysis but not actually a parameter in our numerical method. We explain how we discretize Equations 7 and 10 in Section 4.

Accuracy of approximation. Equation 10 is a linearized approximation to the true evolution, which is valid only when $l$ is close to $k$. Thankfully, the Gaussian term in Equation 9 forces $A$ to be small whenever $l$ is far from $k$, so the errors in $\backslash mathcal\{A\}$ caused by this approximation are exponentially small.

Furthermore, the reflecting boundary condition is based on geometrical optics and is valid only when the wavelength and the packet size (represented by the parameter $s$) are small compared to the boundary curvature. A more accurate handling of high curvature boundaries would require to take wave scattering effects into account, which we leave as future work and note that our simple reflecting boundary condition is sufficient for visual plausibility.

Low-frequency $\backslash mathcal\{A\}$. We prove in Appendix B that the Fourier transform of $\backslash mathcal\{A\}$ is essentially a low-pass-filtered version of the Fourier transform of $\backslash eta$. In other words, we can reconstruct the wave height function $\backslash eta$ using only low-frequency variables. This frequency shift has important consequences in the eventual discretization of the method, because the Nyquist limit prevents us from discretizing the high-frequency $\backslash eta$ function directly without aliasing, while the discretization of $\backslash mathcal\{A\}$ is possible even with coarse grid resolutions. Instead of the grid resolution imposing limits on the resolution of visible water waves, our algorithm will impose limits on the resolution of the wave amplitudes, which are not visualized directly and arguably more difficult to discriminate perceptually.

Phase Shifts. The Gabor transform allows us to discuss the water height in terms of phase as well as local amplitudes. The phase of each wavelet in Equation 7 is determined by $\backslash mathbf\{k\}\; \backslash cdot\; \backslash mathbf\{x\}\; -\; \backslash omega(k)t$. If we translate the entire function in space by replacing $\backslash eta(\backslash mathbf\{x\},\; t)$ with $\backslash eta(\backslash mathbf\{x\}\; -\; \backslash mathbf\{y\},\; t)$, then we gain a phase shift of $-\backslash mathbf\{k\}\; \backslash cdot\; \backslash mathbf\{y\}$ on the right hand side of the equation; a small shift in space can lead to a large phase shift if $k$ is big. Consequently, although the amplitudes $\backslash mathcal\{A\}$ are well-behaved even on coarse grids, the phases are still sensitive to Nyquist's limit and require either low frequency waves (small $k$) or very many samples of $k$ for accurate reconstruction. Alternatively, if we do not require particular interference patterns, then we can avoid aliasing by adding a random initial phase shift to each wave [Cook 1986]:

where $\backslash xi$ is a random number for each wave direction that is fixed throughout time.

Hybrid simulation. Section 3.1 motivates our method as a hybrid between spectrum-based methods and PDE-based methods. This hybrid nature becomes more explicit if we examine the effect of the size of the Gaussian, $s$ on our simulation variable $\backslash mathcal\{A\}$:

When $s$ approaches infinity and the Gaussian becomes a spatial constant, our method transitions to a traditional spectrum-based algorithm. On the other hand, when $s$ goes to zero and the Gaussian becomes a Dirac delta function, our method computes a function similar to the water height, $\backslash eta$. We never set $s$ directly in our discretization, but it is essentially determined by the grid spacing, $\backslash Delta\; x$. Thus, our method resembles a spectrum-based algorithm within a single grid cell, and it resembles a PDE-based discretization as we zoom out.

Relationship to Wave Packets. Although we used the Gabor transform to derive our approach in Section 3.2, we can alternatively derive this method from the water wave packets of Jeschke & Wojtan [2017], as we do in Appendix C. We show that Equation 7 emerges as the limiting behavior of an infinite number of wave packets, spanning all possible positions and wavevectors. From this perspective, Jeschke & Wojtan [2017] introduced one particular Lagrangian discretization of our continuum theory, which samples a small number of individual packets and tracks their propagation through space. Equation 10 represents an alternative Eulerian reference frame, which tracks the changes in the wave packet content at each point in space. We also note that, in Appendix C, the function $\backslash mathcal\{A\}$ appears as a smooth average of all nearby wave packet amplitudes. This gives us further evidence that $\backslash mathcal\{A\}$ is indeed a smooth function that varies slowly over space.

The amplitude $\backslash mathcal\{A\}(\backslash mathbf\{x\},\; \backslash mathbf\{k\},\; t)$ is a function in $4\; +\; 1$ dimensions: two in space, two in wavevector, and one in time. We find it intuitive and computationally convenient to represent the wavevector in polar coordinates $\backslash mathbf\{k\}\; =\; (k\; \backslash cos\; \backslash theta,\; k\; \backslash sin\; \backslash theta)$, where $k$ is the magnitude of $\backslash mathbf\{k\}$, and $\backslash theta$ is the angle made by $\backslash mathbf\{k\}$ and the $x$-axis. We represent $\backslash mathcal\{A\}$ on a four-dimensional grid $[x\_\{\backslash min\},\; x\_\{\backslash max\}]\; \backslash times\; [y\_\{\backslash min\},\; y\_\{\backslash max\}]\; \backslash times\; [0,\; 2\backslash pi]\; \backslash times\; [k\_\{\backslash min\},\; k\_\{\backslash max\}]$, which spans two spatial coordinates $x$ and $y$, the angular coordinate $\backslash theta$, and the wavevector coordinate $k$. We store samples of $\backslash mathcal\{A\}$ on each of the nodes in this 4D grid, indexed by the coordinates $a,\; b,\; c$. We use the notation $\backslash mathcal\{A\}\_\{abc\}$ to represent the discrete amplitude sample of the wave at grid node position $\backslash mathbf\{x\}\_a\; =\; (x\_a,\; y\_a)$, traveling at angle $\backslash theta\_b$, with wavenumber $k\_c$. Figure 2 illustrates the grid used for our discretization.

We can now approximate $\backslash mathcal\{A\}$ with a linear combination of basis functions, in the style of the finite element method:

where $\backslash phi\_a(\backslash mathbf\{x\})$ is a basis function in position, $\backslash theta\_b(\backslash theta)$ is a basis function in angle, $\backslash psi\_c(k)$ is a basis function in wavenumber coordinates, and $C(\backslash mathcal\{A\}\_\{abc\},\; t)$ is a coefficient function weighting various values of

$\backslash mathcal\{A\}\_\{abc\}$. For example, piecewise-linear basis functions in all dimensions would set $C(\backslash mathcal\{A\}\_\{abc\})\; =\; \backslash mathcal\{A\}\_\{abc\}$, giving us

in which the basis functions are simply hat functions that weight the nearest values $\backslash mathcal\{A\}\_\{abc\}$.

The only practical restriction on the basis functions is that our reconstructed function should actually interpolate the discrete samples $\backslash mathcal\{A\}\_\{abc\}$ if we want operations like semi-Lagrangian advection to work as expected; in other words, $\backslash mathcal\{A\}(x\_a,\; \backslash theta\_b,\; k\_c,\; t)\; =\; \backslash mathcal\{A\}\_\{abc\}(t)$. The spatial basis function $\backslash phi\_a(\backslash mathbf\{x\})$ simply interpolates $\backslash mathcal\{A\}\_\{abc\}$ over space, and we use standard piecewise-linear or piecewise-cubic basis functions here. The angular weighting $\backslash vartheta\_b(\backslash theta)$ interpolates amplitudes across different traveling directions, and we use piecewise cubic basis functions here as well.

The wavenumber basis $\backslash psi\_c(k)$ has a special physical interpretation; it symbolizes the spectrum of the waves represented by amplitude $\backslash mathcal\{A\}\_\{abc\}$. If we continue with our wave packet analogy, then $\backslash psi\_c(k)$ describes the shape in frequency-space of the wave packet at position $\backslash mathbf\{x\}$ traveling in direction $\backslash theta$ with representative wavenumber $k\_c$. Thus, a piecewise-constant $\backslash psi\_c(k)$ implies that the wave packet has a flat spectrum, with all wavenumbers similar to $k\_c$ having the exact same amplitude; a piecewise-linear $\backslash psi\_c(k)$ gives the packet a bit more of a localized wave packet-like shape with its peak at $k\_c$; and a Gaussian $\backslash psi\_c(k)$ resembles the typical wave packet derivation. We are free to assign any wave spectrum we wish to this function, and, since we are modeling water waves in this paper, we found it appropriate to use an actual ocean wave spectrum for $\backslash psi\_c(k)$. Our examples set $\backslash psi\_c(k)$ equal to the directional spectrum in Equation 32 of [Horvath 2015], normalized such that $\backslash mathcal\{A\}(x\_a,\; \backslash theta\_b,\; k\_c,\; t)\; =\; \backslash mathcal\{A\}\_\{abc\}(t)$.

In practice, we found that we can use surprisingly coarse grids for this discretization. Our simulations use $\backslash Theta\_\{\backslash mathcal\{A\}\}\; =\; 16$ samples for the wavevector angle, $\backslash theta$; finer discretizations did not increase the simulation quality. We use even fewer samples for discretizing the wavenumber, $k$; we typically only need $K\_\{\backslash mathcal\{A\}\}\; =\; 1$ sample to get the effects we desire, though we experiment with up to $K\_\{\backslash mathcal\{A\}\}\; =\; 4$ samples

in Section 8. As mentioned above, $\backslash mathcal\{A\}$ varies slowly over space, so we do not require much spatial resolution either. In our implementation we allocate $X\_\{\backslash mathcal\{A\}\}\; =\; 4096$ grid cells for each spatial dimension, which defines a grid cell spacing of approximately one meter.

Lastly, the simulations in this paper use real-valued $\backslash mathcal\{A\}$ functions as initial conditions. This function then stays real for all time according to Equation 10, so our implementation does not bother to store complex $\backslash mathcal\{A\}\_\{abc\}$ coefficients.

Once we have discretized all of the wavevectors $\backslash mathbf\{k\}$, Equation 10 becomes a small number of independent scalar advection equations in space—one for each sample of $\backslash mathbf\{k\}$. We numerically integrate each of these equations in parallel using the unconditionally stable semi-Lagrangian advection with slope-limited cubic spatial interpolation proposed by Fedkiw et al. [Fedkiw et al. 2001]:

where $\backslash mathcal\{A\}\_\{bc\}$ is the interpolation of the discretized $\backslash mathcal\{A\}$ with fixed angle $\backslash theta\_b$ and wavenumber $k\_c$, $\backslash hat\{\backslash mathbf\{k\}\}\_b\; =\; (\backslash cos\; \backslash theta\_b,\; \backslash sin\; \backslash theta\_b)$ is the wave direction determined by angle $\backslash theta\_b$, and $\backslash mathbf\{x\}\_a\; =\; (x\_a,\; y\_a)$.

Whenever a semi-Lagrangian ray leaves the simulation domain, we apply the relevant boundary condition from Equation 11 by using a different value on the right hand side of Equation 17. For transmitting boundaries, we assign a procedural $\backslash mathcal\{A\}\_\{\backslash text\{ambient\}\}$ function describing the ambient ocean behavior. For reflecting boundaries, we reflect the semi-Lagrangian ray off the boundary to get a new position $\backslash mathbf\{x\}\_\{\backslash text\{reflect\}\}$ and new direction $\backslash mathbf\{k\}\_\{\backslash text\{reflect\}\}$, and we update $\backslash mathcal\{A\}$ with the reflected function value $\backslash mathcal\{A\}(\backslash mathbf\{x\}\_\{\backslash text\{reflect\}\},\; \backslash mathbf\{k\}\_\{\backslash text\{reflect\}\},\; t)$.

Amplitude Spreading. Amplitudes with different angle $\backslash theta\_b$ or wavenumber $k\_c$ have either different travel direction or speed. Over time these amplitudes separate from each other and a nice initial state can turn into bunch of separated amplitude blobs as depicted on the left part of Figure 3. This is a problem of discretizing Equation 10 with finitely many wavevectors and we deal with it by adding two diffusion terms:

where $\backslash delta$ controls diffusion in the direction of travel (mimicking how wave packets stretch out due to dispersion), and $\backslash gamma$ controls diffusion in the travel angle (mimicking how wave packets spread out tangentially as they radiate from a source). Other diffusive terms are possible, but we found that the above two terms work really well. Figure 3 illustrates the effect of our angular diffusion term.

We apply operator splitting to solve this advection-diffusion equation numerically: the advection is solved by semi-Lagrangian advection as described above, and the diffusion is discretized with second order finite differencing in space and forward Euler in time. In our examples, we set the diffusion parameters relative to the advection speed $\backslash omega\text{'}$, the spatial resolution $\backslash Delta\; x$, the wavenumber resolution $\backslash Delta\; k$ and the angular resolution $\backslash Delta\; \backslash theta$: $\backslash gamma\; =\; 0.025\; \backslash omega\text{'}(k)\; \backslash Delta\; \backslash theta^2\; /\; \backslash Delta\; x$ and $\backslash delta\; =\; 10^\{-2\}\; \backslash Delta\; x^2\; \backslash Delta\; k^2\; |\backslash omega\text{'}\text{'}(k)|$.

To calculate the actual water height we numerically evaluate the integral in Equation 12. We evaluate the wavenumber in polar coordinates

and plug in the discretized form of $\backslash mathcal\{A\}$ to get:

Now we can see more explicitly how the basis function $\backslash psi\_c(k)$ relates to the spectrum associated with amplitude $\backslash mathcal\{A\}\_\{abc\}$; Equation 21 treats $\backslash psi\_c(k)$ as the amplitude function for each wavelength $k$ in the heightfield evaluation.

To use the above formulas, we first pre-compute the function $\backslash bar\{\backslash Psi\}\_c(p,\; t)$ at the beginning of each time step by evaluating it at a number of discrete sample points $p\_i$ and store it in a 1D texture, which we call a “profile buffer”. This pre-computation allows us to turn this costly integral evaluation into a simple texture lookup for the rest of the time step. After the pre-computation, we can evaluate the water height at any point in space by approximating the one-dimensional integral in Equation 20 with a summation over several sampled angles. To create our examples, we sum over $K\_\backslash eta\; =\; 400$ wavenumbers for the pre-computation in Equation 21, and we store the resulting function $\backslash bar\{\backslash Psi\}\_c$ in a 1D texture of size $N\; =\; 4096$. We sum over $\backslash Theta\_\backslash eta\; =\; 120$ angles to evaluate $\backslash eta$ at each point in space. Note that the number of directions and wavenumbers for discretizing $\backslash mathcal\{A\}$ compared to $\backslash eta$ are largely independent. We chose these parameters empirically based on the quality of our results, and we list them in Table 1.

Our examples also extend the water height field to mimic trochoidal Gerstner waves [Tessendorf 2004b] by computing horizontal displacements in addition to the vertical ones in the profile buffer $\backslash bar\{\backslash Psi\}\_c(p,\; t)$. We then include these horizontal displacements in the final summation of $\backslash eta$.

The pre-computation of the profile buffer in Equation 21 is largely responsible for the performance of our method. For comparison, the evaluating $\backslash eta$ at all points in space using the naive 2D integral in Equation 19 costs $O(K\_\backslash eta\; \backslash Theta\_\backslash eta\; X\_\backslash mathcal\{A\}^2)$ operations, where $X\_\backslash mathcal\{A\}^2$ is the total number of spatial locations where $\backslash eta$ is computed. Our speedup using the profile buffer reduces this computation by two orders of magnitude to $O(K\_\backslash eta\; +\; \backslash Theta\_\backslash eta\; X\_\backslash mathcal\{A\}^2)$. Distributing this computation over $G$ GPU cores reduces the cost further to $O(K\_\backslash eta\; +\; \backslash Theta\_\backslash eta\; X\_\backslash mathcal\{A\}^2/G)$. In practice, the introduction of the profile buffer in our implementation raised the frame rate from 1.8 to 275 – a speed-up factor of 233 for evaluating $\backslash eta$.

1: function TIMESTEP($t$)
2:   AdvectionStep($t$)
3:   WavevectorDiffusion($t$)
4:   $\backslash bar\{\backslash Psi\}\; \backslash leftarrow$ PrecomputeProfileBuffers($t$)
5: end function

6: function WATERHEIGHT($\backslash mathbf\{x\},\; t$)
7:   $\backslash eta\; \backslash leftarrow\; 0$
8:   for $b\; \backslash leftarrow\; 1,\; \backslash Theta\_\backslash eta$ do
9:     $\backslash theta\_b\; \backslash leftarrow\; \backslash frac\{2\backslash pi\}\{\backslash Theta\_\backslash eta\}\; b$
10:    $\backslash hat\{\backslash mathbf\{k\}\}\; \backslash leftarrow\; (\backslash cos\; \backslash theta\_b,\; \backslash sin\; \backslash theta\_b)$
11:    $\backslash mathbf\{p\}\; \backslash leftarrow\; \backslash hat\{\backslash mathbf\{k\}\}\; \backslash cdot\; \backslash mathbf\{x\}\; +\; \backslash text\{rand\}(b)$
12:    for $c\; \backslash leftarrow\; 1,\; K\_\backslash eta$ do
13:      $\backslash eta\; \backslash leftarrow\; \backslash eta\; +\; \backslash mathcal\{A\}(\backslash mathbf\{x\},\; k\_c\; \backslash hat\{\backslash mathbf\{k\}\})\; \backslash cdot\; \backslash bar\{\backslash Psi\}\_c(\backslash mathbf\{p\},\; t)$
14:    end for
15:  end for
16: end function

In nature, wave energy is lost due to viscosity, splashes and so on. In order to reproduce this behaviour, we add the same dissipation term as Jeschke et al. [2017] to Equation 18:

The first term models dissipation due to water viscosity $\backslash nu$ ($10^\{-6\}\; \backslash text\{m/s\}$), which affects mostly waves with small wavelengths. The second term accounts for surface contamination such as oil, dirt or algae and it affects larger waves as well.

The method of Jeschke & Wojtan [2015] computes steady-state wave motions by pre-computing the paths of multiple wavefronts, storing the phase functions on an adaptive triangle mesh, and proposing a phase interpolation scheme to recover the final waves at runtime. We can produce similar results with our simulator by computing

wave motion until the amplitudes reach a steady state, and then storing the $\backslash mathcal\{A\}$ functions as static textures. We can then interpolate these $\backslash mathcal\{A\}$ textures at any point in space and use them to efficiently compute wave heights at run-time. This pre-computation of $\backslash mathcal\{A\}$ leads to a considerable speed-up over running the full simulation ($4.7\backslash times$ speedup from $\backslash approx\; 60\backslash text\{fps\}$ to $\backslash approx\; 280\backslash text\{fps\}$ in this example). Our method is arguably more efficient than Lagrangian wavefront tracking (though the methods have very different numerical parameters and the errors behave differently), and the run-time behavior is more convenient and load-balanced for graphics hardware than an adaptive triangle mesh.

Figure 4 shows an example scene that was pre-computed using this technique. To make this result, we modeled wind effects by adding a source term $S(\backslash mathbf\{x\},\; k,\; t)$ to Equation 18. Local wind effects like this are difficult for wavefront tracking methods which require the explicit creation of coherent wavefronts as initial conditions, but they are easy for our Eulerian method. However, unlike [Jeschke and Wojtan 2015] our simulator does not yet handle shallow water effects due to a spatially varying dispersion relation $\backslash omega$, and our jittered sampling strategy makes it difficult to precisely control the wavefront phases.

We implemented some elementary coupling between our fluid simulator and rigid bodies. We do this coupling by adding forces to the rigid body from the fluid, and adding waves to the fluid from the rigid body each time step.

To make the waves influence the rigid bodies, we compute a buoyancy force by approximating the shape of the submerged volume of the rigid body as a cylinder with volume $V\; =\; \backslash pi\; r^2(\backslash eta\; -\; r\; -\; h)$, where $r$ is the radius of the body and $h$ is the height of its center of mass. The buoyancy force is then $-\backslash rho\; V\; \backslash mathbf\{g\}$, where $\backslash rho$ is the density of the fluid and $\backslash mathbf\{g\}$ is the gravity vector. We integrate this force using backward Euler integration to ensure numerical stability.

To add waves to the rigid body simulation, we calculate the change in energy of the rigid body caused by the water $\backslash Delta\; E\_\{RB\}$, where the rigid body's energy $E\_\{RB\}$ is equal to $m\backslash mathbf\{v\}^2/2\; +\; mgh$, with body mass $m$, velocity magnitude $v$, and gravity magnitude $g$. Linear wave theory tells us that water wave energy in deep water is proportional to the squared amplitude:

and we use our localized amplitudes to confine the wave energy to a region of space:

To approximately conserve energy, we must set $\backslash Delta\; E\_\{\backslash text\{fluid\}\}\; =\; -\backslash Delta\; E\_\{\backslash text\{RB\}\}$. Assuming this change in energy is distributed equally among all wavevectors, we can then compute the new amplitudes $\backslash mathcal\{A\}$ at the location of the rigid body $\backslash mathbf\{x\}\_\{\backslash text\{RB\}\}$:

where the product $K\_\{\backslash mathcal\{A\}\}\; \backslash Theta\_\{\backslash mathcal\{A\}\}$ is the number of discrete wavevectors in our discretization.

To summarize, our rigid body coupling is executed each time step by adding buoyancy forces to each body, calculating the change in energy of the body, and then converting that change in energy into waves at the location of the rigid body. Figure 5 uses this technique in an example with numerous floating boxes. This heuristic coupling strategy works well for displaying basic solid-fluid interactions, and we leave a more careful treatment for future work.

The basis function $\backslash psi(k)$ controls the wave spectrum which is visualized in the final results. The spectrum can either be determined by physics or tuned by hand to create more stylized results. Figure 6 shows how changing this function affects the visualized wave heights. Regardless of the chosen spectrum, the waves will travel at the correct phase speed due to the dispersion relation in Equation 21.

Instead of using physical equations to compute the amplitude textures $\backslash mathcal\{A\}$, we can explicitly create or modify them with procedural functions, or with a novel amplitude-painting interface. In our video, we show that we can procedurally create reflection effects by artificially amplifying waves traveling perpendicular to an obstacle's surface, with an amplification factor based on distance to the obstacle. We can also use a painting interface, as shown in Figure 7.

We note that hand-tuning these amplitudes will be physically incorrect in the sense that the group speeds are set to zero instead of determined by the dispersion relation. However, the phase speeds are still physically correct due to the dispersion relation in Equation 21. The result is that the waves themselves travel at the correct speed, but the wave groups do not. We found this indirect physical inaccuracy a bit more difficult to perceive, and so we believe that such an amplitude override technique might be a useful artistic tool.