Abstract
Simulating large amounts of water in real time is often achieved using heightfield methods. Allowing the water to interact with rigid bodies is essential for applications such as games, but traditional heightfield interaction methods concentrate on water-to-body effects by letting water flow through the bodies. We instead take an approach where the bodies block water. Our earlier method is improved in several ways, taking steps toward a single method to create both water-to-body and body-to-water effects. The new method is also visually compared to a traditional method by Thürey et al. A drawback of our method is that it has some grid aliasing artifacts that appear especially when the method is used for floating bodies. However, our method is demonstrated to work together with the Thürey method, which allows us to get the best of both worlds to simulate both floating and blocking bodies in a single scene. The method runs in real time for large areas of water even with a very limited GPU budget.
1. Introduction
Water is a common and important element in nature. It is also frequently encountered in games and other 3D virtual environments. A lot of attention has been devoted to rendering water realistically, but interaction with it is typically limited in games. Bodies of water are often simply modeled as static planes with possibly some procedural waves that have no effect on gameplay and therefore no interactivity. Interaction is sometimes possible with small amounts of particle-based liquids, but they are also more commonly used only for visual effect.
Interactive systems are at the heart of gaming. Off-the-shelf rigid body physics solvers have revolutionized 3D environments in part because they provide so many natural interaction possibilities. Similar interaction with large amounts of water is currently not possible in the real-time 3D realm. A more versatile interaction of rigid body physics and large-scale bodies of water could enrich virtual worlds tremendously and even enable completely new game genres. In the 2D domain, this has already been achieved by games such as Where’s My Water [1] and Sprinkle [2] using particles. Some first steps have also been taken in 3D by the pioneering game From Dust [3].
Performance is one of the main reasons why large-scale water areas are still static in most game worlds. Large-scale fully 3D water is still out of reach, but the continuing increases in the available parallel GPU computing power are making some simpler methods fast enough for consideration.
Both Lagrangian (particle-based) fluid solvers, such as smoothed particle hydrodynamics (SPH) [4], and fully 3D Eulerian (grid-based) solvers have been extended to handle interaction with solids in very advanced and realistic ways, with current research focusing more on the pathological cases. In the particle-based fluids, even the real-time results are very impressive for small amounts of fluids [5]. Sadly, neither of these approaches are fast enough for real-time modeling of large bodies of water, such as rivers and lakes [6].
Heighfield-based methods simplify the situation by giving water a single depth value at each 2D coordinate point. This allows them to be much faster than particles or fully 3D methods for large water areas, but a lot of surface detail is lost [7]. From a gaming perspective, an even bigger problem is that interaction with rigid bodies is currently very limited in the heightfield context.
To make interaction possible, water needs to affect the bodies via buoyancy and drag. On the other hand, the bodies also need to affect the water, because they block the flow of water, causing waves, and may even make a river take a completely different route.
None of the previous heightfield-based methods solve both interaction directions in a satisfactory way. The traditional solution in real-time methods is to let the water be mostly unaffected by the object, that is, allowing water and bodies to overlap (e.g., [8–14]). In these methods, buoyancy and drag are easy to apply based on the submerged part of the object. Some surface waves may then be generated to create an illusion of the objects affecting the water. However, in these models, water flows through the bodies.
In real world, the objects would naturally push down the water surface. Unfortunately, this breaks the heighfield assumption as soon as water should flow on top of the body. This can be worked around with some limiting assumptions on the blocker shapes or a multilayer model. However, there is hardly any research in this direction, and the few previous methods either have had to let most of the water pass through the bodies [11] or have not been able to include water-to-body effects [15]. The latter method also has problems with rendering.
In this paper, we improve and generalize the method of [15]. This enhances the effect of water on the objects, solves the rendering issues, takes steps toward more general blocker geometry, and replaces some ad hoc solutions by a more physically based approach. The approach and its remaining problems are also analyzed in much more detail than in [15].
Unfortunately, the problem of spatial aliasing from the grid that has plagued the Eulerian 3D fluid-body interaction methods is also present in our method. These problems are mostly negligible for heavy objects but become clearly visible if the method is used to handle floating bodies. As an alternative solution for applications where the aliasing is found to be too disturbing, we propose combining the two interaction approaches. Objects are classified as floaters or blockers. This allows both body-to-water and water-to-body effects in a single scene, as in Figure 1.
All presented methods run faster than real-time on relatively large terrains using our GPU-parallel implementation.
The rest of the paper is structured as follows. Section 2 introduces some related work. Our simulation method, including body-to-water effects, is described in detail in Section 3. Section 4 discusses the problems of adding water-to-body effects in it and suggests a combined method. Section 5 evaluates the results and Section 6 concludes.
2. Related Work
Fluid simulation has a long tradition in engineering. However, most methods from this literature strive for realism, for example, bridge building, and are much too slow for games. In computer graphics, water simulation started to become popular in the 1990s with, for example, the work of Foster and Metaxas [16], but the main application has been in special effects for movies, with processing times for single frames often measured in minutes. We concentrate on methods that are suitable for real-time purposes. For a good introduction on the more realistic approaches, see Bridson’s book [7] for the Eulerian approach and the recent review of the SPH literature for the Lagrangian perspective [17].
2.1. Fluid Simulation Methods
Lagrangian (particle-based) methods, such as SPH, are currently popular in the research community [17]. Because the particles are fully 3D, particles are needed to reach a given spatial resolution, which is not yet fast enough for modeling rivers, lakes, and other large-scale features in real time [18]. On the other hand, only the areas containing water need to be simulated, and thus performance is limited by the amount of water, not the total area or volume. Creating a renderable surface also takes a long time compared to heighfield methods [18]. Despite this, particles have been successfully used in 3D game engines for small-scale effects such as leaky pipes [19] and in 2D games even on mobile platforms [1, 2]. Particle-based fluid solvers also have the advantage of having many sophisticated and robust methods for solid-fluid interaction [5, 20–22], though most of them are not designed for real-time applications.
Eulerian methods track the fluid quantities in a fixed grid. While fully 3D methods are mostly too slow for games, there are also faster approaches. Tall cells only model a thin layer near the water surface in 3D, and the underlying noninteresting part is treated as a single, tall cell [23, 24]. This approach yields very high-quality results even in real time but is not yet fast enough for general usage in games, because they have multiple other systems requiring computational resources [6].
The Lattice Boltzmann method has also been used to simulate water [25]. To our knowledge, it is not fast enough for the task of simulating large-scale water in real time due to small time steps. The Lagrangian and Eulerian approaches can also be combined. For example, in movie effect production, FLIP is currently a very popular method, which combines aspects of both methods [26]. As a fully 3D method, it is also too slow for most real-time games.
In open seas without a need for dynamic flow over rough and partly dry terrain, wave particles [12] and FFT-based methods [27] can be used [18]. They can provide relatively realistic-looking results and interaction with rigid bodies, for example, ship wakes [10].
The most relevant methods for our discussion are based on heighfields. Only allowing a single height value for the water surface removes many interesting surface effects, such as breaking waves and splashes, but has the great advantage of making the simulation for a given spatial resolution. Furthermore, these methods have no need for an expensive linear system solver unlike the Eulerian simulations. Thus heighfield methods allow for the simulation to cover vast areas in real time. Many of the missing surface phenomena can be created using fast procedural methods and only near the viewer, because they are mostly aesthetical and not essential for gameplay.
The pipe method is a very simple and fast water simulation method [8, 28]. It has been extended for multilayer situations [29] and used for purposes such as erosion simulation [30, 31]. Many GPU implementations have also been presented, demonstrating its speed on modern hardware [6, 30, 32].
Shallow water equations provide a more physically based approach to heightfield water simulation [7]. They are slightly more complicated than the pipe method but still fast enough for games. They have also been extended in various ways and implemented on the GPU [11, 13, 33]. However, there is some evidence that the visual difference to the pipe method is quite small when modeling large-scale water bodies in a gaming context [34].
2.2. Interaction with Solids
Simulating the interaction of fluids and solids has been the focus of much research in the fully 3D fluid simulation literature. Lagrangian methods typically model the solids also as particles and are able to produce impressive results even in real time [5]. However, our Eulerian method does not have much in common with these methods.
In the Eulerian context, Foster and Metaxas voxelized the boundary conditions to the simulation grid [16]. Rigid bodies were later handled by Takahashi et al. [35] by rasterizing the rigid body velocities to the grid and setting these as velocity border conditions to the fluid. Many methods run the fluid and rigid body solvers alternatingly, but the quality can be improved by solving both problems simultaneously, as in [36].
The approach of rasterizing border conditions to the grid has problems with boundaries that are not aligned to the grid [37]. A large body of research has addressed these problems since then. A common approach based on the immersed boundary method [38] is to average solid properties in each cell (or dual cell) of the grid and use these as weights in the solver [37, 39]. Some authors also started to model the solids directly in a Lagrangian framework instead of rasterizing them to the grid [36]. A different approach, mostly unsuited for real-time applications, is to give up the grid regularity either completely or at least near the solids [23, 40].
The fully 3D methods have lately been applied to real-time situations, but only on rather small grids. Many of the methods can just as well be applied on the 2D grids, which would be large enough for many games, but this requires the solids to also reside in a 2D domain. In a heightfield situation, the water and the rigid bodies still reside in a three-dimensional domain, but the fluid simulation grid is two-dimensional. Most methods, such as that of Batty et al. [37], are based on modifying the pressure projection step, which does not even exist in heightfield simulations. Because of these reasons, the sophisticated 3D fluid-body interaction methods cannot be directly applied when the underlying simulation is based on heightfields. Thus the heightfield methods have traditionally taken a completely different, and much simpler, approach to the problem.
The heightfield water-body interaction research has mostly concentrated on handling floating objects [8–10, 14]. In these methods, the bodies are superimposed on the water, which makes buoyancy and drag almost trivial to calculate. The effect of the bodies on the water is handled by, for example, slightly pushing down the surface at the bodies, which creates believable waves and ripples for floating objects. The problem of cell aliasing is also present in these methods, and the iWave method solves it similarly to, for example, [37] by estimating an opacity value for the blocker coverage in each cell [10]. Some researchers even calculate interaction per mesh triangle, subdividing the triangles if necessary [12, 13].
However, the previously mentioned heightfield methods are unvariably only interested in the interaction of the solids with the water surface, completely ignoring the body of the water. Approaching water waves are passed through the bodies almost unaffected, which makes sense for small and floating objects that only move up and down with the waves. On the other hand, large or heavy bodies are not handled in a sensible way. This problem is most evident if a body is lying on dry ground and a wave meets it (see Figure 8(a), or the video in the Supplementary Material available online at http://dx.doi.org/10.1155/2014/580154), or when one tries to build a dam out of some heavy bodies. A method that is also able to handle submerged objects was introduced in [11]. It is still mostly interested in surface effects and does not block flow. A method that allows large bodies to block flow by forcing water to flow above, below, or around the body was introduced in [15] but had several shortcomings that we try to address here.
3. Simulation Method
This section describes our simulation method, which is a variant of the classic pipe method of [8] with a layer structure similar to that introduced in [29]. This work is a continuation to [15].
We use a top-down approach, first introducing the general structure of our model in Sections 3.1 and 3.2. A detailed description of each phase follows in Sections 3.3–3.6. Finally, we discuss the parameters and stability in Section 3.7.
3.1. Model Structure
Water behavior is simulated on a uniform 2D grid on the -plane. The distance between neighboring cell centers in both directions is . The model consists of cell-sized columns of terrain, water, blocker, and air stacked on top of each other. A column occupies the vertical interval . The depth of a column is the vertical length of the interval it occupies. A column with zero depth is considered not to exist. Blocker columns are nonpenetrable to water and are used to represent the moving rigid bodies, which are rasterized to the grid resolution.
The model is generalizable to allow any number of water and blocker columns stacked on top of each other, but in practice we have found it enough to have a single terrain and a single blocker column plus up to two water columns (one above the blocker, one below), with an implicit air column above each water column. A more general structure would greatly increase the memory bandwidth requirements. The bottom column is always a heightfield terrain that occupies the interval . Instead of storing the top and bottom coordinates of water columns, a single height (in meters) is stored for each. When needed, the coordinates can be calculated based on the terrain and body locations. The water surface is a heightfield that is constructed by finding the top of the topmost nonempty water column at each cell.
The blocker columns represent the rigid bodies and are found using the rasterization hardware as will be described in Section 3.3. A cell is free, if there is no blocker column. Free cells only have a single water column and a single air column. The water columns are combined by removing the air column between them. A water column is free if there is no body above it. Figure 2 illustrates the structure of our model.
Water is allowed to flow between neighboring cells in the grid. The four von Neumann neighbors are used, but speed could be traded for quality by using the Moore neighborhood. If both cells have a blocker, we assume the blockers come from the same body. This is reasonable, because it is extremely rare that different bodies are located in neighboring cells but are not overlapping. Overlapping bodies would be treated as one anyway, because we only have a single blocker column.
The previous assumption allows us to only store two flows for each pair of neighboring cells: “upper” and “lower.” If both cells are free, upper contains the flow, while lower is 0. If exactly one of the cells is free, the two stored flows are free top and free bottom. If neither of the cells is free, upper contains the flow between the top layers and lower between the bottom layers.
We define as the set of neighbors of water column . The flow between columns and is (naturally only one value is stored). If , it is called outflow from cell and inflow otherwise.
The external pressure at the top of each water column is also included in the model. External pressure is the result of a rigid body over the column and is zero for free water columns. Air pressure is ignored as negligible. The pressure inside a water column is assumed to behave according to the basic hydrostatic equation, , where is the pressure at the top of the column, is water density, and is acceleration by gravity. We always give pressure as , that is, depth (in meters) of a water column that would cause the pressure. The pressure caused by a rigid body above a water column is modeled as explained in Section 3.6.
Since in this method water pushes the surface down, buoyancy and other forces affecting the bodies cannot be calculated directly in the usual way. Due to the finite grid size, rendering the surface directly would also cause disturbing dents near the objects as illustrated in Figure 3. To solve these problems, an effective surface is constructed to approximate where the surface would be if there were no blockers. The effective surface is another heightfield with value equal to the actual surface plus the external pressure at the surface (in meters). Figure 3 illustrates the concept.
3.2. Algorithm Phases
The following phases are executed to advance the simulation for a single time step of :(1)rigid body update (by an external physics engine),(2)blocker update (apply effects of the moved bodies to the water, Section 3.3),(3)flow update (Section 3.4),(4)depth update (Section 3.5),(5)pressure estimation (Section 3.6).
Each step is described in more detail below, but let us first give an overview of the steps. Our method is based on alternating rigid body and water update steps. The rigid body update is a standard simulation step of an off-the-shelf physics engine, Bullet [41] with buoyancy, and other fluid forces applied to the bodies. After that, the water overlapping with the moved bodies is removed, as described in Section 3.3.
Somewhat analogously to the simplifying assumption of hydrostatic pressure always used in the pipe method, we also assume that the rigid bodies are static when solving for water flow and depth. The solid interaction in flow and depth update steps (Sections 3.4 and 3.5) can thus be seen as generalizations of how terrain is typically handled in the corresponding steps of the pipe method, only with several layers and thus the possibility of water flowing both above and under the blocking body. This departs from the more physical approach of the 3D methods, where solid velocities are taken into account, but we have found the results to be quite acceptable.
The final step, pressure estimation (Section 3.6), is needed for visualization and estimating buoyancy and drag.
The method is designed to be completely parallel on the GPU with no random write access. Each phase uses double buffering, only reading the previous state. This ensures that the fast texture memory can be used.
For clarity, no time indices are used for the variables. All formulas are given in an assignment form instead.
3.3. Blocker Update
The blockers are rigid bodies that have possibly moved in phase 1 of the algorithm, which will affect the water.
For each cell, the rasterization finds the bottom and top height coordinate of any rigid bodies overlapping it at the center. This is implemented by rasterizing the triangle meshes using an orthogonal view from above and below. The depth buffer can then be employed to find the extremes, which makes the process very fast even for complicated bodies. See Figure 4 for a visualization of the concept.
In the rasterization step, multiple bodies located on top of each other are combined to a single blocker column. Bodies completely above water are ignored. Otherwise, water between an object even high above the water and a submerged object would be removed.
The next step is to ensure that the bodies and water do not overlap. The blocker columns before and after the time step are handled as being completely unrelated. Note that the vertical extent in a cell blocked by the same body can change radically even during a single time step as the body moves and rotates.
In each cell, the procedure is implemented as follows. The possible blocker column that existed before the time step is first removed. Water is added so that the surface is moved to where the effective surface was. Then the possible new blocker column is added, removing any overlapping water. The addition divides the water column to (possibly empty) parts below and above the blocker. The external pressure in the column below the body also needs to be updated so that the effective surface is not changed (see pressure approximation in Section 3.6).
During the whole process, the depth of artificially added and removed water is kept track of. Their difference could be either positive (water was added) or negative (water was removed). To conserve volume, the same amount needs to be removed from or added to the vicinity. This is achieved by averaging the amount of missing or extra water with the four neighbors on each time step. Whenever a free cell has a value indicating missing or extra water, a portion of the difference is added or removed. This is what mostly causes waves around a body dropped into water. Our approach is very similar to that used in [11].
An example of this phase can be found in Figure 5. For more details, see [15], which has a similar implementation.
3.4. Flow Update
In addition to removing water from inside the bodies, it is also essential to prevent any flow from entering the bodies. To achieve this, we let the cross section of the pipes connecting columns to vary, similarly to [29], but taking blocker geometry into account.
The flow between two water columns and that are neighbors goes through an interface that is represented by the red arrows in Figure 2. The bottom of the interface . For the top, we first interpolate the surface height at the interface by taking the average of the tops of the water columns (dashed lines in the figure). We also find the bottom of the lowest blocker above the interacting water columns, which sometimes limits the interface (e.g., lower flow between cells and in Figure 2). Putting these together, the top of the interface is . The cross section of the interface is .
Flow is caused by a pressure gradient. In column , the pressure relative to the bottom of the interface is , with calculated similarly. The pressure difference is , analogously to the original pipe method of [8].
The flow can now be updated using an explicit Euler step as in the standard pipe method: where is a friction coefficient.
Finally, the flow from each column is restricted to prevent the depth from becoming negative. This is done by calculating a scaling factor , where is the net outflow. Each outflow from is then multiplied by .
3.5. Depth Update
The depth of each column can now be updated based on the flows. The flow through each interface in a single time step is limited to the depth of the interval, . This step is not physically based but takes care that intervals that have become almost or completely blocked do not let water through due to previously accumulated flow. It is therefore a key element of making the bodies block flow. For example, water cannot flow on top of a body before the surface is high enough. This is roughly analogous to how terrain is handled in the related methods (e.g., the reflective boundary conditions caused by high terrain in [13]). After this limiting, the update assignment is
Finally, some water columns are under a blocker. If this update would cause a water column to overfill so that it would overlap the blocker above, only water equal to the free portion is allowed to flow. Experimentally, we found out that stability is greatly increased by also decreasing the incoming velocities to such a cell, which can be thought of as a crude model of a no-slip border condition with the blocker. If is the depth of the overlap and is the total incoming depth (not net incoming depth) during a time step, we multiply all inflows to the cell by . This approach also removes the artifacts caused by the need to slow down underflow in [15].
3.6. Approximating Pressure
The effective surface is needed for buoyancy and rendering but is not directly available in places where a body pushes the surface down. Figure 3 is referred to as an example throughout this section to help understand the situation.
It would be possible to track the pressure under a body using physical simulation and calculate the effective surface that way. Instead, we use a simpler approach, and interpolate the effective surface from around the bodies. The external pressure is zero in the free columns, which allows the surface to be used directly as the effective surface (cells , , and through in the example). These are the border conditions of the interpolation.
In practice, the interpolation is achieved by applying a blur filter on the effective surface of the previous time step and storing the external pressure calculated from that. The calculated values are only used in the nonfree columns. In Figure 3, the results are applied to cells –. As a few time steps go by, the effective surface in these cells gradually rises from the surface, approaching the dashed line in the figure.
However, this process is complicated by the bottom water columns (in and ). Their external pressure could be calculated directly as the difference to the surface. Now imagine the water level at changing slightly so that the top water column in it dries up. This would suddenly push the effective surface below the body, which would result in visually disturbing flickering and unstable buoyancy. We work around this by keeping track of the pressure also in the bottom water column of such cells and letting the external pressure be the maximum of the directly calculated value and the iteratively solved value.
Formally, let be a nonfree column and the cell it is located in. The current external pressure , stored in each such column , is updated during a time step: where is the set of neighboring cells of and controls the spreading speed of pressure.
The maximum in (3) affects cells such as, for example, , where is the top of the top water column. While the body in is underwater, the external pressure keeps the effective surface at , but now, after dries up, the external pressure stays. As a result, the effective surface changes very smoothly, which is essential to keep the visual appearance continuous in time.
For a visual comparison of the surface and the effective surface, see Figure 6 and the accompanying video.
(a)
(b)
3.7. Parameters and Stability
As a variant of the pipe method, our solver also uses explicit time integration. This sets a limit to the time steps that can be used before disturbing instabilities are experienced.
The proposed method has a few parameters that affect the stability and other properties of the simulation: , , , , and . The stability properties are hard to analyze exactly due to the complexity of the blocking, but the idea of the CFL condition [7] also carries over to our method: to make the method stable, a short time step is needed to compensate for large wave speeds (where wave speed is measured in simulation cells per second).
A longer time step makes the simulation use less computational resources. While the time step might affect the simulation in some subtle ways, a previous user study did not find any noticeable effect when varying the time step inside a stable range in a similar simulation without the rigid body interaction [34]. Therefore, after fixing the values of the other parameters, the time step should be chosen as large as possible without making the simulation unstable.
The friction parameter, , describes how much of the old flow is conserved between time steps and should therefore be a function of the time step length. A smaller value makes the water appear more viscous but also limits the velocities accumulated during the simulation, therefore allowing longer time steps. The simulation scale is controlled by . A smaller (a finer grid resolution) needs to be compensated with a smaller time step or smaller velocities.
A larger causes water displaced by the moving bodies to be added back faster, causing bigger waves to emanate from a blocking body in the water. This has the side effect of causing spikes and large flow velocities if the value is too large. Finally, a larger makes the gaps near objects disappear faster but also makes the effective surface look more twitchy, since it reacts faster to blocker geometry changes.
Typical values for the parameters in our simulations are , ms, m, , and . With these values, the stability problems have been minimal.
4. Water-to-Body Effects
The method described thus far does not implement any effects caused by the water to the bodies. These effects can be divided into two classes: drag and buoyancy. By drag we mean that water and object velocity tend to equalize. This is especially important in the case of floating bodies, such as logs, that need to follow the flow. Both of these effects are easy to achieve if the water surface is allowed to go through the floating object. In our model, the actual surface is pushed down by the object, but the effective surface can be used in its stead. Our approach here is fairly standard but is shortly described for the sake of completeness.
For the drag calculation, we additionally need a velocity in each cell . Exactly as in [30], the component can be estimated by calculating the average flow in the direction divided by depth of the column in question and correspondingly. Only the top column is used since drag is mostly relevant for floating bodies.
Drag is based on iterating through all cells overlapped by the body. For each cell , we first calculate , which is the the vertical extent of the body below the effective surface. We also calculate the local velocity of the rigid body in the center of the cell at the middle of the underwater part. The drag force , where is a drag coefficient that can be set empirically according to the needs of the application. In our scenes, . The force is applied at the middle of the underwater part during each time step.
Buoyancy is equal to the weight of the water displaced by the body, again applied in each cell overlapped by the body. We reuse from above and for each cell apply at the center of the cell.
One problem in this approach is that the drag and buoyancy forces are calculated in only the grid locations, using relatively long time steps, which often causes excess rotation. This can be fixed by simply damping both the linear and rotational velocities of the bodies using an empirically determined coefficient based on the proportion of the body that is underwater. This works relatively well in practice but is not physically based. However, our implementation of the method presented in [11] also needed to resort to similar tricks.
4.1. Aliasing Problems
The suggested blocking method (and that of [15]) is binary: each cell is either blocked or nonblocked. As explained in the related work section, the problem is a common theme in the Eulerian fluid simulation literature. In our case, where the quality goal is much less ambitious, this problem is usually negligible. However, a very small change in the body position can cause relatively large changes in the rasterized version, when a body suddenly overlaps a cell center. This spatial aliasing can cause overly large effects to the water that become noticeable if the body is only moving slowly.
The blocking method is thus best suited for situations where the bodies mainly move independently of the water. If a body stops completely, the physics updates can be disabled and the aliasing problem disappears. On the other hand, applying even an otherwise insignificant buoyancy force to a heavy object keeps the body alive, which causes small vibrations in most rigid body solvers and makes the problem visible. How disturbing one finds this phenomenon depends on the scene and parameters.
Unfortunately, one of the worst cases is an object floating along a strong current. When a cell is suddenly blocked, the stable flow to the cell is denied, causing the source cell to be quickly overfilled. This small movement causes a comparatively large wave to emanate from that cell to the surroundings. This looks unnatural, since there is almost no relative movement between the body and the water, meaning no waves should be created. The problem is illustrated in Figure 7.
(a)
(b)
(a)
(b)
Several directions for a solution were examined in our research process. Simply limiting the maximum velocity or making the grid denser makes the effect smaller but does not solve the problem without making the method too slow. The best simple workaround we found is limiting the amount the surface in any cell can change in a single time step to some constant, but it only makes the problem somewhat more apparent.
A better solution could be based on the fractional cell method, such as in [37] or [10]. However, they are not trivially applicable, as trying to take the fully 3D nature of the problem into account is problematic with our underlying 2D fluid solver and trying to address the full complexity of the 3D situation would quickly become too slow. A strictly 2D solution such as in [10] also does not work as such in our context, because our method completely relies on the vertical intervals to block the flow. In reality, those intervals vary in complex ways inside the cell, and a single opacity value is not enough to describe them, which makes the problem difficult to solve.
4.2. Combined Solution
Because a single solution that would implement both water-to-body and body-to-water effects in all kinds of scenes in a satisfactory way is not yet available, we propose an alternative implementation that combines the strengths of the two approaches. Bodies are categorized according to their density to floaters and blockers. Floating bodies are advected with the flow and do not really need to block water. They can therefore be treated with the Thürey et al. method [11]. On the other hand, bodies that are denser than water do not float and can be treated with our method without the previously mentioned aliasing problems becoming an issue.
The Thürey et al. method only affects water in a single phase that changes the depth of the topmost water column. This can be implemented between steps 1 and 2 in our method without any adverse interaction with our method. Water-to-body effects in our examples are implemented as suggested above, using the drag and buoyancy forces, but also, for example, the triangle subdivision scheme of [13] could be used.
5. Results and Discussion
To evaluate and illustrate our results, we have created three interactive test scenes. The first two demonstrate blocking capabilities. In the first one, the user drops heavy cubes and two types of cars on a smooth hill with a large water source on top (Figures 8 and 9). The second demonstrates controlling river flow by dropping blocking obstacles (Figures 11, 12, 13, and 14). The last one combines floating with blocking. It has logs floating down a river with a car crossing the river, hitting some of the logs on the way (Figures 1 and 10). The scenarios use either a or a simulation grid with a 1-meter simulation resolution. The time step is 20 ms.
We have elected to benchmark our method by comparing it to our implementation of the interaction method of [11], which we find to be the most advanced of the real-time heightfield methods, especially when it comes to handling heavy and large bodies that can be submerged. Other heightfield methods such as [9, 13, 14] would look very similar to the Thürey et al. method in our test cases with heavy objects, since none of these prevent water from passing through the bodies. For the floating objects, we do not claim that our method currently achieves a quality similar to any of the mentioned methods. A comparison to the significantly more advanced methods based on a 3D grid or particles is not included, since these cannot be directly applied to a heightfield-based simulation. None of the heightfield-based methods are currently even close to being competitive with them in quality, but, on the other hand, there is a large performance benefit when using heightfields.
The first two scenarios are handled well by our method. Water flows around the blocking rigid bodies and under a truck. Visual comparisons between our method and the one in [11] have been provided in Figure 8, Figures 11 through 13, and the accompanying video. The differences should be very apparent, since the Thürey et al. method [11] lets water mostly flow through the bodies, only creating limited interaction waves that are in many cases dwarfed by the existing water movement.
The method of [15] also fares well in these cases. The main difference is in our pressure tracking, which allows for much nicer object/water borders than those achieved in [15].
In the third scenario, the logs can in principle be implemented as blocking bodies using our method. They float and are advected by the flow, but the spatial aliasing problem as discussed in Section 4.1 makes the results unconvincing. This is one of the main limitations of our method. The advection is also somewhat limited, since the velocity under the body is usually zero. This could possibly be solved in the future by interpolating the velocity field for the effective surface similarly to pressure.
When the third scenario is implemented with the combined solution described in Section 4.2, the logs are classified as floaters and do not block water. The car is a blocker and is not affected by the water. We find this scenario to work much better with the combined solution than when implemented with either of the methods alone.
The stability of our method was tested empirically by varying the time step. It was found that, even with high-speed objects, the simulation is stable as long as the underlying pipe method itself is stable, that is, the rigid body interaction is not what limits the selection of the time step. The only direct effect our interaction method has on fluid velocities is decreasing them. The blocker update sometimes creates large height gradients when a cell is left empty, but we did not find any noticeable stability problems in practice even in those cases. In the combined solution, the first problem appearing when increasing the time step was the floating objects (whether they used our method or that of [11]) becoming unstable due to the torque caused by the buoyancy forces. Finally, when it comes to the aliasing problems, the time step does not play much of a role.
Our method and its predecessor [15] are the only heightfield-based interaction methods where flow is really blocked by the bodies. The other methods are concentrated on scenarios such as floating boats, where the body-to-water effect is mostly intended for visual, rather than interactive, purposes. However, we think that controlling water masses with rigid bodies will open up many kinds of new interaction possibilities in games and should be given more consideration. On the other hand, we also demonstrated that these two kinds of effects can coexist in a single scene.
The method of [15] also allows for this kind of interaction but is less grounded in physical simulation than our blocking solution. Most importantly, we take external pressure into account similarly to [29] and take steps toward a more general approach that would allow the method to be extended to any number of layers and more complex geometry. For example, new physical situations can now be simulated, such as the equalization of two columns of water that are connected by a tunnel, because the external pressure in the tunnel is included in the model.
A crucial addition compared to [15] in our method is the procedure for estimating the effective surface. The result is roughly the same surface that would exist in, for example, the method of [11] and can therefore be used for calculating water-to-body effects. It also has the necessary property of reacting smoothly to all kinds of geometry changes, which is necessary for retaining visual quality.
Thürey et al. [11] use a similar idea to our blocker handling step as the primary method for handling object-to-water effects. This has the problem that water is moved blindly to the surroundings, causing water to leak through a blocking body. In our approach, most of the work is done in a later step when calculating the new flow, which limits flow into the bodies.
In all examples, the water simulation and body-to-water effects are implemented on the GPU. The model is designed for implementation without random write access, which allows using fast GPU texture memory without any need for thread synchronization or atomic operations. Our implementation packs the simulation state into four grid-sized textures with four 32-bit floating point channels each (and their duplicates for double buffering).
On a laptop with an NVIDIA Quadro 1000 M, a simulation step takes about 5 ms of the GPU time on a grid and 2 ms on a grid. For the 1 m resolution grid in our examples, a 20 ms time step is stable enough. This makes the method feasible for real-time applications with large bodies of water, even with the limited GPU budget that this kind of system would be allowed in a game.
Solving the aliasing problems could result in a single method usable for all kinds of scenarios. A priority for future work is to either generalize the 2D approach of [10] to our more complicated case or modify and simplify a 3D approach such as the one in [37] to become compatible with our underlying 2D simulation and fast enough for large grids. Our method is also admittedly not rigorously based on physics. Momentum conservation is not considered, the solid-fluid velocity equality constraint is not enforced, and the depth update still has some ad hoc steps to improve stability. Finding a way to apply a completely physically based solution without losing too much speed is another important topic for future work. Our combined solution could also be improved by dynamically changing between the methods depending on the situation. For example, if a floating piece of ice gets stuck, it could be converted to a blocking body, allowing for a dam to be formed.
Our implementation concentrates on the water behavior and interaction and is visually fairly simple. Implementing adaptive tessellation, spray, foam, bubbles, breaking waves, and other effects using, for example, a particle system would naturally improve the visual quality of the results considerably. It could also be possible to extend the method to create, for example, foam based on the actual interaction, instead of the typical texture-based procedural methods used in current games (similarly to [42]).
6. Conclusions
We presented a GPU-parallel real-time water-solid interaction method for heightfield water simulation. Our method extends and generalizes that presented in [15]. It is based on bodies pushing down the water surface but includes a procedure for estimating where the surface would be without the body. This allows richer rigid body interaction in heightfield water simulation than what could be achieved using any of the previous methods.
A limitation of our method is that some aliasing artifacts are encountered especially when the method is used to handle floating objects. Despite trying various approaches, the artifacts could not be removed. To circumvent the problem, combining different methods in the same scene was proposed and demonstrated.
One of the most important remaining problems is to find a solution for the spatial aliasing problems inherent in the presented flow blocking method. This could allow a single method to be used for all kinds of bodies. Further physical realism should also be strived for, but it is important to do this without sacrificing performance. Various visual enhancements are also possible, including foam, and are integrating a particle system to create splashes.
The presented method is fast enough for real-time games on current systems, even with the need to dedicate most of the computational resources for various other game functionalities. The method can be combined with an off-the-shelf rigid body solver, enabling new kinds of player interaction with water.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was partly funded by KAUTE foundation. The author would also like to thank the editor and anonymous reviewers for their useful comments that helped to improve the quality of this paper.
Supplementary Materials
The supplementary material contains a video that illustrates the water-object interaction method described in the article. The method is also compared to a previous one.