Abstract
A variety of numerical techniques have been explored to solve the shallow water equations in real-time water simulations for computer graphics applications. However, determining the stability of a numerical algorithm is a complex and involved task when a coupled set of nonlinear partial differential equations need to be solved. This paper proposes a novel and simple technique to compare the relative empirical stability of finite difference (or any grid-based scheme) algorithms by solving the inviscid Burgers’ equation to analyse their respective breaking times. To exemplify the method to evaluate numerical stability, a range of finite difference schemes is considered. The technique is effective at evaluating the relative stability of the considered schemes and demonstrates that the conservative schemes have superior stability.
1. Introduction
Fluid dynamics is widely used in computer game technology to simulate a range of phenomena including water, smoke, and soft bodies. A variety of numerical techniques have been explored for simulating fluids in real-time such as position-based dynamics (e.g., [1, 2]), the finite-volume method (e.g., [3–6]), and the finite-element method [7]. There are two broad modelling approaches to describing fluids: grid-based Eulerian or particle-based Lagrangian methods [8, 9]. The Lagrangian approach follows an individual fluid parcel as it moves through space and time, whereas the Eulerian approach describes the fluid motion from specific fixed points in space as time passes. Both approaches have enjoyed widespread use in game technology. Particle-based approaches are advantageous when considering arbitrary boundaries, fluid mixing [10], and interactions with rigid bodies [11]. Grid-based approaches are widely used in computer graphics applications and can attain higher numerical accuracy since dealing with spatial derivatives is easier to accommodate on a fixed grid. However, in contrast to particle-based methods, grid-based approaches have difficulty ensuring the conservation of mass and can be computationally slower. Furthermore, in real-time water simulations in games, grid-based methods perform much better at tracking smooth water surfaces [9].
Smooth Particle Hydrodynamics (SPH) is a popular particle-based method for simulating fluids in game technology, computer animation, virtual reality, and the movie industry. For example, incompressible SPH is a promising numerical scheme for large-scale and large-deformation simulations used in interactive fluid flow simulations [12]. Two of the challenges faced by SPH is unstable solid boundary handling and numerical dissipation, both of which inhibit stable and realistic fluid evolution. Using a position-based dynamics framework integrated into the SPH solver [2, 13], these issues can be overcome [14]. SPH is also useful for simulating viscoelastic materials such as gels, gelatin, and mucus by applying a velocity correction to limit the fluid deformation, making it attractive to soft-body simulation as well. This produces visually accurate results but is offset by computational performance costs; however, simulations can be accelerated using GPUs [15]. Recently, neural networks have been used to accelerate particle-based fluid simulations that run significantly faster than a GPU Position-Based Fluid Solver whilst preserving visual quality [16].
The Lattice-Boltzmann method (LBM) is a grid-based method that indirectly solves the fluid dynamical equations by solving the Boltzmann equation that describes the underlying particle distribution function [17–19]. The Boltzmann equation is easier to numerically solve in contrast to the classical fluid equations, and LBM can be run efficiently on massively parallel architectures. LBM is useful for solving complex, coupled, and multiphase flows and can accommodate complex boundaries with relative ease [17]. However, tracking and preserving small-scale features, such as a fluid drop or splash, is a challenging problem. Recent work has looked at a novel grid-particle method for reconstructing distribution functions of interface grids and coupling the reconstruction method with the LBM and Volume of Fluid method to track the free surface [20]. The method enhances the accuracy of the reconstruction and helps preserve the fluid surface detail.
The shallow water equations (SWEs) are a simplified set of fluid flow equations that are widely used in computer graphics applications. The SWEs are well suited to game technology applications due a number of practical reasons: in contrast to the full 3D Navier-Stokes equations, the relative simplicity of the SWEs lead to a significant performance advantage; the solutions produce a full velocity field which is useful for providing plausible interactivity and solid-fluid coupling, and the rendering of the simulation results can be undertaken with common, fast rendering techniques [8, 21]. In contrast, full 3D simulations require surface tracking, volumetric rendering techniques, or mesh generation algorithms for visualisation. Although the SWEs can describe complex nonlinear fluid phenomena such as vortices, it is unable to handle breaking waves or splashing phenomena [8, 21].
A variety of numerical methods exist for solving the SWEs. Finite difference (FD) techniques are particularly popular in real-time applications due to the relative simplicity of their implementation and their ability to capture complex phenomena such as shocks. In particular, a range of finite difference techniques has been applied to the conservative form of the shallow water equations including the Lax-Friedrichs (LF) scheme, the Richtmyer two-step Lax-Wendroff (LW) method, and the McCormack (MC) method [22, 23]. Variations of these popular algorithms have also been investigated such as the corrected LF (CF); composite schemes, e.g., LWLF, MCLF, and CFLC [24, 25]; filtered LW, MC, and CF schemes; and the Picard Integral Formulation of the Weighted Essentially Non-Oscillatory (WENO) scheme [26]. Many of the numerical schemes, such as the LW, MC, and CF methods, suffer from oscillatory phenomenon (high-frequency jitter) as a result of instability. As a consequence, attempts are made to mitigate instabilities from occurring. The filtered schemes tackle the stability issues by using an oscillation smoothing method, for example, Liang et al. [27] and Hsieh et al. [28] for stabilizing their tsunami simulations using the MC scheme. Similarly, Ransom and Younis [29] used a total-variation diminishing limiter with the MC method to help avoid oscillatory behaviour, hence stabilizing the scheme.
Recent studies have proposed a generalized finite difference-split coefficient matrix method along with the flux limiter technique to eliminate potential numerical oscillations that could lead to instability [30]. Li et al. presented a fifth-order weighted essentially nonoscillatory scheme for simulating dam break flows in a finite difference framework [31] that produces smaller truncation errors and provides the same accuracy order and stability as contemporary WENO schemes. To accommodate more complex systems, such as flat-bottom geometry, a new discretization of the source term of the SWEs is suggested by Prieto-Arranz et al. [32]. In the approach, a Smooth Particle Arbitrary Lagrangian-Eulerian formulation based on Riemann solvers is used to solve the SWEs where stability is achieved using a posteriori MOOD paradigm. Benchmarking demonstrates that the MOOD limiting procedure is able to prevent artificial oscillations occurring in the neighbourhood of discontinuities and shocks. Wu et al. presented a high-order entropy stable discontinuous Galerkin method for solving the SWEs on curved triangular meshes that preserves a semidiscrete entropy inequality and remains well-balanced for continuous bathymetry profiles [33]. Such an approach is advantageous as it can accommodate complex geometries through unstructured meshes; it is simple to parallelize and is able to take advantage of acceleration techniques using GPUs.
Determining the stability of a numerical algorithm is a complex task. For linear cases, the von Neumann stability analysis [22] can be applied analytically but is intractable in nontrivial scenarios including complex nonlinear problems, where a set of coupled partial differential equations have to be solved. In nontrivial scenarios, elaborate numerical tests can be used (such as circular dam breaks and flows around a bump) to give an in-depth characterisation of the numerical algorithms performance (e.g., see Parna et al. [26]), but this approach can be time-consuming, complex to set up, and difficult to analyse.
This paper proposes a novel and simple technique to compare the relative empirical stability of finite difference algorithms (or any grid-based scheme) by solving the inviscid Burgers’ equation to analyse their respective breaking times. A similar technique has been used to determine suitable plasma fluid solvers in astrophysics [34]. The proposed technique provides a quick and easy way of determining the stability of a proposed algorithm prior to more thorough stability tests tailored to the specific system to be solved and its application. The proposed method is in keeping with the method of A- and L-stability where numerical methods for stiff problems involving ordinary differential equations are analysed by applying them to the test equation for , . In the fluid context, the inviscid Burgers’ equation emulates the form of partial differential equations (PDEs) in which there is an advective term such as in Euler’s equations of fluid dynamics and the SWEs.
The inviscid Burgers’ equation (or the nonlinear wave equation) is given by
It describes a nonlinear wave propagating in the positive -direction with a speed proportional to its amplitude, . Therefore, larger wave amplitudes propagate faster than smaller components: this causes the wave to break due to nonlinear steepening, like the phenomenon of breaker ocean waves, yielding an advective instability. Obtaining the characteristics for the first-order PDE reveals a solution that may become triple-valued after a period of time—this is wave breaking. The onset of wave breaking is characterised by the solution having a vertical tangent at the leading wave edge. The timescale for this to first occur is the breaking time. In the numerical context, multivalued solutions are not possible and the numerical schemes become unstable. The ensuing nonlinear instability manifests itself as a series of oscillating spikes that form in the wake of the leading wave edge. The relative empirical stability of two numerical methods can be assessed by analysing their respective breaking times. Due to the nature of the instability, where a discontinuity forms due the solution becoming triple-valued, the method also gives an indication of how well numerical schemes deal with discontinuities such as shocks and boundary interactions. To exemplify the method to evaluate numerical stability, a range of FD schemes will be considered.
This paper is structured as follows. The Methodology introduces and derives the FD algorithms that will be used to exemplify the stability analysis technique. It outlines the initial and boundary conditions for the stability evaluation and calculates the corresponding analytical solution, including the breaking time and location, for the system. The Methodology concludes by defining a quantitative metric to evaluate the numerical stability of the numerical schemes. Following this, the Results and Discussion evaluates the numerical stability of the FD algorithms using the quantitative metric. In this paper, a comparison is made of the relative stability of FD numerical methods, but other grid-based numerical schemes could be similarly analysed.
2. Methodology
The finite difference methods are popular numerical schemes for solving partial differential equations (PDEs). The utility and efficacy of these numerical methods depend upon a number of criteria including stability, accuracy, and ease of implementation.
The LW method transforms a continuous problem into a discrete one by replacing spatial and temporal derivatives with second-order accurate finite difference approximations, thus reducing the problem to an iterative algebraic exercise [22]. The LW method is reliable, stable, and accurate. However, depending on the complexity of the individual PDE or if the system is governed by a coupled set of PDEs, implementing LW can be very complicated.
Finite Difference Time Domain (FDTD) is an alternative numerical method [35–37]. It entails recasting the PDE as a series of ordinary differential equations (ODEs) with respect to time by discretizing the spatial domain and replacing the spatial derivatives with finite difference approximations. The ODEs are then solved numerically via an ODE solver such as the Runge-Kutta or leap-frog methods. Its main advantage lies in its ease of implementation.
2.1. Lax-Wendroff Method
To obtain a finite difference solution of (1), the domain described by the PDE is defined in terms of a rectilinear grid, its edges parallel to the and axes. We define discrete coordinates so that an arbitrary grid point is given by such that and , where , are integers and , are the grid spacings in the , coordinates, respectively. Writing , the finite difference formulae are obtained from the Taylor expansion of about in the neighbourhood of , while holding fixed [22],
This is second-order accurate in time. The LW method contains an inherent diffusion term, , that helps dampen the instability. For the case of the linear wave equation, is replaced by and by , where is a constant wave speed. Therefore,
Using the finite difference approximations: This yields the Lax-Wendroff algorithm for the linear wave equation, where is the mesh ratio. This algorithm is second-order accurate in space and time. One of the advantages of the finite difference method is that it is very easily extended to the solution of nonlinear equations since in all but special cases many of the methods and proofs are invariant [22].
For the nonlinear equation (1), the finite difference formula (3) is used where is replaced by , using (1), and is replaced by the following,
Substituting into the finite difference expression (3) yields
The Lax-Wendroff algorithm for the nonlinear wave equation becomes where the finite difference approximations (5) and (6) have been used.
Numerical algorithms for solving partial differential equations are only useful if they are convergent and stable. A finite difference algorithm is deemed convergent if the difference between the theoretical solutions of the differential and difference equations at a fixed coordinate tends to zero when the number of grid points in a fixed size numerical domain is increased: and . An algorithm is considered stable when the difference between the theoretical and numerical solutions of the difference equation remains bounded as tends to infinity. The von Neumann stability analysis for the linear wave equation algorithm states that the scheme is stable provided , which appropriately coincides with the Courant-Friedrichs-Lewy (C.F.L) condition for the convergence of the algorithm. This result is used as a guide for the stability and convergence of the nonlinear expression; ergo, [22].
An alternative numerical algorithm is possible if the second righthand term of (10) is differenced directly: substituting (10) into (3) gives where
This yields an alternative algorithm
Note that this version uses next nearest neighbours. Additional finite difference algorithms can be derived by considering the conservative form of the nonlinear wave equation,
Following the previous derivation and using (3), the term is replaced with
Substituting (19) into (3) and using yields where . The Lax-Wendroff algorithm for the conservative form of the nonlinear wave equation is
Using (18) in (3) instead of (19) yields
Differencing the derivative in the final term directly produces an alternative conservative algorithm:
Therefore,
2.2. Finite Difference Time Domain Method
The FDTD method is as follows. Using the predefined discrete coordinates and writing , the spatial derivative in the nonlinear equation was replaced by the central difference formula (5), yielding
This is a first-order ODE with respect to time where is the distance between successive spatial grid points. Solving using a 4th-order Runge-Kutta method gives where is the time step, is the number of time steps, and is the integration time such that . By considering the conservative form of the nonlinear wave equation, another FDTD algorithm can be derived. Substituting for the central difference formula yields
The corresponding Runge-Kutta algorithm is, therefore,
Both the FDTD algorithms described here are second-order accurate in space and fourth-order accurate in time. For stability, the algorithms require that , where is the maximum expected phase velocity. This ensures that the solution cannot vary significantly over one spatial increment during one temporal step [36].
2.3. Initial and Boundary Conditions
To compare the FD algorithms, they were used to solve the nonlinear wave equation, the inviscid Burgers’ equation,
For the initial conditions, the wave amplitude was perturbed from a uniform equilibrium in the shape of a Gaussian waveform
, where is the amplitude of the perturbation, defines its width, and is the centre of the computational domain (see plot of Figure 1). Each algorithm was computed for a prescribed length of time for the specified initial condition. So that the LW and FDTD methods could be compared, the values of and had to be chosen carefully. Setting requires that implying . For the simulations, the boundary conditions were for all . For the simulations discussed here, the parameters listed in Table 1 were used, without loss of generality.
2.4. Analytical Considerations of the Inviscid Burgers’ Equation
To evaluate the relative empirical stability of the FD algorithms, it is instructive to evaluate the breaking time and position analytically for the prescribed initial and boundary conditions. The inviscid Burgers’ equation is given by
Performing the change of variables and yields where . For the initial conditions, a Gaussian wave packet is prescribed:
A solution to this initial value problem can be sought using the method of characteristics. The characteristics are where is a constant on a given characteristic. The implicit solution for is
One can rewrite the solution as an explicit function for ,
A plot of versus for various evaluated using (41) is plotted in Figure 2. The solution shows the wave amplitude as a function of position for various values of up to, including and beyond, the breaking time.
The solution (Equation (40)) may become triple-valued after a period of time—this is wave breaking. The onset of wave breaking is characterised by the solution having a vertical tangent at the leading wave edge. The breaking time is the time at which this first occurs. Recall the characteristics, where . Defining , one can find a particular value of such that . Using , the breaking time is simply
Given the breaking time , one can calculate the corresponding breaking location, , by solving the characteristics for with and ,
Therefore, the breaking location
For the particular problem described in this manuscript, parameterized by the values in Table 1, the breaking time and location are and , respectively (illustrated by the solid curve in Figure 2). In comparison, for the numerical solutions, the wave will break earlier () and at a premature location , which is to be expected since the numerical grid resolution will be defeated long before the leading wave edge approaches an infinite gradient.
2.5. Evaluating the Numerical Stability
To compare and quantify the development of the instability, we define the fractional change of the numerical solution relative to the analytical solution, where is the analytical solution (Equation (40)) evaluated at and , the breaking time and breaking location of algorithm , respectively, and is the numerical solution of algorithm , evaluated at and . The earlier the instability evolves for a particular algorithm, the smaller the value of and the poorer the algorithm performs. The bigger the value of , the poorer the stability of the algorithm. Note that the definition of is consistent with the quantity, , defined by Mitchell and Griffiths [22], to determine the stability of a finite difference numerical algorithm.
To determine the breaking time, , and breaking location, , for algorithm , for each time step, , the maximum value of the wave amplitude, , is determined for . Using this, the absolute magnitude of the difference between the maximum amplitude and the amplitude of the initial perturbation, , is calculated: . In the first instance that this quantity is greater than or equal to a threshold value, , then the corresponding values of and are the breaking location, , and breaking time, , of algorithm , respectively. In this paper, the threshold value, , is set to be 1% of the amplitude of the initial perturbation, .
3. Results and Discussion
In the discussion that follows, LW1, LW2, LW3, LW4, FDTD1, and FDTD2 denote numerical solutions of (1) obtained from the formulae (13), (16), (21), (24), (31), and (33), respectively. Figure 1 shows the evolution of the nonlinear wave equation solved using FDTD1. The plot encapsulates qualitatively the typical behaviour of the nonlinear wave equation as a function of position and time. The initial, perturbed wave-form propagates in the positive -direction with a speed proportional to its amplitude. Therefore, small amplitude wave components propagate slower than their larger counterparts, leading to wave steepening. In general, as fast moving wave elements overtake those moving slower, the wave solution becomes multivalued for a single value of . Within the numerical context discussed here, multivalued solutions are not possible and the numerical schemes become unstable. The ensuing nonlinear instability manifests itself as a series of oscillating spikes that form in the wake of the leading wave edge. As the leading wave edge tries to overtake further wave elements, more spikes are formed in its wake. Figure 3 exhibits the initial development of the advective instability using FDTD1 and FDTD2 as example cases. The plot shows a close-up of the wave crest for time steps (cross), (plus sign), and (point), illustrating the formation of the initial spike that disrupts the solution.
(a)
(b)
To quantify the development of the instability, we calculate the fractional change of the numerical solution relative to the analytical solution, . For the algorithms LW1, LW2, LW3, LW4, FDTD1, and FDTD2 considered here, , 0.0412, 0.0273, 0.0279, 0.0473, and 0.017, respectively (see Table 2), demonstrating that in order of stability from best to worst, we have FDTD2, LW3, LW4, LW2, LW1, and FDTD1. The nonconservative algorithms appear to have inferior stability in comparison to the conservative algorithms. Of the nonconservative algorithms, the FDTD1 solution has the least stability, while the two nonconservative LW methods (LW1 and LW2) are virtually identical. Of the conservative algorithms, the FDTD2 solution has the most stability, and very little separates the two LW solutions (LW3 and LW4). Although the FDTD2 algorithm appears to break before the LW algorithms (as indicated by its value), its comparison to the analytical solution for the same position and time indicates that it ultimately has superior stability. Note that can also be viewed as a measure of accuracy, since instability leads to inaccuracy and so is a dependent concept.
4. Conclusion
This paper proposes a novel and simple technique to compare the relative empirical stability of finite difference algorithms (or any grid-based scheme) by solving the inviscid Burgers’ equation to analyse their respective breaking times. The proposed technique provides a quick and easy way of determining the stability of a proposed algorithm prior to more thorough stability tests tailored to the specific system to be solved. The relative empirical stability of two numerical methods can be assessed by analysing their respective breaking times. Due to the nature of the instability, where a discontinuity forms due the solution becoming triple-valued, the method also gives an indication of how well numerical schemes deal with discontinuities such as shocks and boundary interactions. The proposed method is analogous to the A-stability method but tailored to PDEs in which the advective derivative is key.
The technique is effective at determining the relative stability of grid-based numerical algorithms. It demonstrates that the conservative schemes behave very similarly to one another as do the nonconservative schemes and that the stability of the conservative algorithms is marginally better than the nonconservative algorithms. In order of most to least stable, we have the conservative FDTD algorithm (FDTD2), the conservative Lax-Wendroff algorithms (LW3 then LW4), the nonconservative Lax-Wendroff algorithms (LW2 then LW1), and the nonconservative FDTD algorithm (FDTD1). The LW algorithms contain an inherent diffusive term that can help to dampen the instability, quenching the wave steepening for a while. In contrast, the FDTD algorithms lack a diffusive term to help counteract the inevitable instability. In spite of this, the FDTD algorithms do well to sustain coherent behaviour. In principle, one could add a stabilizing diffusive term:
This is Burgers’ equation where is the viscosity coefficient. If the diffusive term quenches the wave steepening, a stable solitary wave structure persists. Alternatively, the time-fractional inviscid Burgers’ equation is worth investigating: the inclusion of a fractional time derivative may stave off the breaking time for longer, potentially improving stability [38]. Fractional-order differential equations (FDEs) have gained importance and popularity recently in their application to physical systems due to their nonlocal properties. This means that the next evolutionary state of the system depends on its historical states not just its current state. Recent investigations have analysed the time-fractional Navier-Stokes equation [39] and the SWEs [40] for bespoke systems. Many numerical schemes have been proposed for solving FDEs [41, 42] as well as for those applicable to fluid simulation, e.g., the fractional Burgers equation [43], the fractional diffusion equation [44] (relevant to the incompressible Navier-Stokes equation), and the fractional parabolic differential equations [45] (relevant to vorticity-stream function formulation of fluids).
The efficacy of a numerical algorithm depends upon a number of criteria including stability, accuracy, and ease of implementation. LW methods are, by definition, second-order accurate in space and time. In general, the choice of the ODE solver in FDTD dictates the temporal order of accuracy: here, a fourth-order Runge-Kutta method was used. However, other techniques can be used such as the leap-frog method. In this respect, FDTD has superior versatility since greater temporal accuracy can be achieved. Additionally, FDTD can be easily modified to have a higher spatial accuracy by simply using higher order finite difference approximations of the spatial derivatives used. Extending the LW method in such a way is nontrivial. In comparison to the LW method, FDTD is very easy to implement, particularly when a system of coupled differential equations is considered. For a system of equations, the LW algorithm can become hard to obtain as the Jacobian of the system becomes more involved. With this in mind, it would seem reasonable to opt for FDTD for certain problems, especially when a system of differential equations is being considered, such as for the numerical solution of the shallow water equations.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
CRS gratefully acknowledges support from the Division of Games Technology and Mathematics at Abertay University.