A continuous/discontinuous Galerkin solution of the shallow water equations with dynamic viscosity, high-order wetting and drying, and implicit time integration

A continuous/discontinuous Galerkin solution of the shallow water equations with dynamic viscosity, high-order wetting and drying, and implicit time integration

Simone Marras smarras@stanford.edu Michal A. Kopera Emil M. Constantinescu Jenny Suckale Francis X. Giraldo Stanford University, Dept. of Geophysics, Stanford, CA 94305, U.S.A. Naval Postgraduate School, Dept. of Applied Mathematics, Monterey, CA 93943, U.S.A. Argonne National Laboratory, Mathematics and Computer Science, Argonne, IL 60439, U.S.A. The University of Chicago, Computation Institute, Chicago, IL 60637, U.S.A.

The high-order numerical solution of the non-linear shallow water equations (and of hyperbolic systems in general) is susceptible to unphysical Gibbs oscillations that form in the proximity of strong gradients. The solution to this problem is still an active field of research as no general cure has been found yet. In this paper, we tackle this issue by presenting a dynamically adaptive viscosity based on a residual-based sub-grid scale model that has the following properties: it removes the spurious oscillations in the proximity of strong wave fronts while preserving the overall accuracy and sharpness of the solution. This is possible because of the residual-based definition of the dynamic diffusion coefficient. For coarse grids, it prevents energy from building up at small wave-numbers. The model has no tunable parameter.

Our interest in the shallow water equations is tied to the simulation of coastal inundation, where a careful handling of the transition between dry and wet surfaces is particularly challenging for high-order Galerkin approximations. In this paper, we extend to a unified continuous/discontinuous Galerkin (CG/DG) framework a very simple, yet effective wetting and drying algorithm originally designed for DG [Xing, Zhang, Shu (2010) Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations Adv. Water Res. 33:1476-1493]. We show its effectiveness for problems in one and two dimensions on domains of increasing characteristic lengths varying from centimeters to kilometers.

Finally, to overcome the time-step restriction incurred by the high-order Galerkin approximation, we advance the equations forward in time via a three stage, second order explicit-first-stage, singly diagonally implicit Runge-Kutta (ESDIRK) time integration scheme. Via ESDIRK, we are able to preserve numerical stability for an advective CFL number 10 times larger than its explicit counterpart.


1 Introduction

The shallow water equations (Saint-Venant Saint-Venant1871 ()) are a common approximation to the –dimensional Navier-Stokes equations to model incompressible, free surface flows. Due to the ability of high-order Galerkin methods to keep dissipation and dispersion errors low (Ainsworth et al. ainsworthEtAl2006 ()) and their flexibility with arbitrary geometries and hp-adaptivity, these methods are proving their mettle for solving the shallow water equations (SW) in the modeling of non-linear waves in different geophysical flows (Ma ma1993 (), Iskandarani iskandaraniEtAl1995 (), Taylor et al. taylorTribbia1997 (), Giraldo giraldo2001 (), Giraldo et al. giraldoHesthavenWarburton2002 ()), Dawson and Aizinger dawsonAizinger2005 (), Kubatko et al. kubatkoEtAl2006 (), Nair et al. nairThomasLoft2007 (), Giraldo and Restelli giraldoRestelli2010 (), Xing et al. xingZhangShu2010 (), Kärnä et al. karnaEtAl2011 (), Hendricks et al. hendricksKoperaGiraldo2015 (), Hood karolineTHESIS2016 ()). Although it is typically assumed that high-order Galerkin methods are not strictly necessary, they do offer many advantages over their low-order counterparts. Examples include their ability to resolve fine scale structures and to do so with fewer degrees of freedom, as well as their strong scaling properties on massively parallel computers (Müller et al. mullerEtAl2016 (), Abdi et al. abdiEtAl2016 (), Gandhem et al. gandhamEtAl2015 ()). High-order methods are often attributed with some disadvantages as well. For example, they are constrained by small time-steps. To overcome this restriction, we follow Kärnä et al. karnaEtAl2011 () and implement an implicit Runge-Kutta scheme based on Giraldo et al. giraldoEtAl2013 (). Furthermore, the numerical approximation of non-linear hyperbolic systems via high-order methods is susceptible to unphysical Gibbs oscillations that form in the proximity of strong gradients such as those of sharp wave fronts (e.g. bores). Filters such as Vandeven’s and Boyd’s vandeven1991 (); boyd1996 () are the most common tools to handle this problem, as well as artificial diffusion of some sort. However, we noticed that filtering may not be sufficient as the flow strengthens and the wave sharpness intensifies. We have recently shown this issue with the solution of the non-linear Burgers’ equation by high order spectral elements in (marrasEtAl2015b, , §5) and will show how a proper dynamic viscosity does indeed improve on filters for the shallow water system as well. To preserve numerical stability without compromising the overall quality of the solution, Pham Van et al. phanVanEtAl2014 () and Rakowsky et al. rakowskyEtAl2013 () utilized a Lilly-Smagorinsky model lilly1962 (); smagorinsky1963 (). To account for sub-grid scale effects, viscosity was also utilized in the DG model described by Gourgue et al. gourgueEtal2009 () to improve their inviscid simulations. Recently, Pasquetti et al. pasquettiGuermondPopoICOSAHOM2014 () stabilized the high-order spectral element solution of the one-dimensional Saint-Venant equations via the entropy viscosity method. Michoski et al. michoskiEtAl2016 () compared artificial viscosity, limiters, and filters for the (modal) DG solution of SW, concluding that a dynamically adaptive diffusion may be the most effective means of regularization at higher orders. The diffusion model that we present in this work stems from the sub-grid scale (SGS) model that was proposed by Nazarov and Hoffman nazarovHoffman2013 () to stabilize the linear finite element solution of compressible flows with shock waves. This was later applied to stabilize high-order continuous and discontinuous Galerkin (CG/DG) in the context of stratified, low Mach number atmospheric flows by Marras et al. in marrasEtAl2015b (). This approach is based on the ideas of scales splitting in large eddy simulation, where the scales of physical importance are split into grid resolvable and unresolvable. The unresolved scales are then parameterized via the SGS model at hand (Dyn-SGS). Dyn-SGS is applicable to any numerical method and is not limited to Galerkin approximations; its implementation in existing numerical models is straightforward. Stabilization becomes more important as the grid resolution is decreased. The metric that we utilize to evaluate its importance is the spectrum of kinetic energy across wave numbers. To understand the impact of diffusion with respect to the energetic behavior of the numerical solution, we compare the spectra obtained for the stabilized continuous and discontinuous Galerkin solutions against their inviscid counterparts. As viscosity is accounted for, the solution improves at coarser resolutions for both CG and DG. In the case of DG, however, dissipation is already built-in by the flux computation across elements. For this reason, the amount of artificial diffusion that DG requires for stabilization is smaller than the one required by CG. This property of DG has been exploited to build implicit Large Eddy Simulation (ILES) models (e.g., Uranga et al. urangaEtAl2011 ()). ILES relies on the numerical dissipation of the discretization method to model the effects on the unresolved, sub-grid scale eddies rather than using an explicit sub-grid scale (SGS) model, under the requirement of high-resolution rider2006 (); drikakisEtAl2007 (). The MPDATA method of Smolarkiewicz smolarkiewicz1983 (); smolarkiewicz1984 () is an example of such an ILES model used for modeling geophysical flows.

Finally, there is a known difficulty in including wetting and drying algorithms (typically designed for low-order methods) while preserving high-order accuracy. The application of wetting/drying with discontinuous Galerkin using low-order Lagrange polynomials can be found in, e.g., Bunya et al. bunyaEtAl2009 (), Vater et al. vaterEtAl2015 (), Kärnä et al. karnaEtAl2011 (), or Gourgue et al . citegourgueEtal2009, and using Bernstein polynomials up to order three in Beisiegel and Behrens beisiegelBehrens2015 (). The positivity preserving limiter of Xing et al. xingZhangShu2010 () (Xing-Zhang-Shu limiter from now on) was designed for high-order discontinuous Galerkin to solve this problem in particular. Because it is mass conservative and preserves the global high order accuracy of the solution, we implemented it in our model and extended it to continuous Galerkin as well. This limiter is guaranteed to work if and only if the mean water height is positive. Its designers recommend using the total variation bounding (TVB) limiter of Shu shu1987 (); cockburnShu1989 () to make the solution positive definite. Because the TVB limiter requires a bound on the derivative of the solution, which we do not know a priori, we rely on the artificial dissipation introduced above instead.

Looking at the larger picture of things, the ever increasing interest in the high-order accuracy of inundation models stems from a long list of coastal disasters in the last 10 years alone. From the 2004 tsunami in the Indian Ocean that caused 280,000 deaths, to the 2011 Tohoku tsunami, followed by Superstorm Sandy in 2012, the devastating 2013 typhoons in the Philippines and India, and the 2014 hurricane Odile; the list seems to be getting longer rapidly as the global climate patterns are changing. To study the impact of similar events in the future, coastal planners around the world are more and more relying on numerical tools such as the one that we present in this paper.

The rest of the paper is organized as follows. The shallow water equations are defined in §2. Their numerical solution is described in §3. The inundation model is described in §5 and is followed by the derivation of the SGS model in §4. Numerical results are reported in §6. The conclusions are drawn in §7.

2 Equations

Let be a fixed domain of space dimension with boundary and Cartesian coordinates in 1D and in 2D and let identify time. Let always identify the direction of gravity. The absolute free surface level and bathymetry are identified by the symbols and, respectively, , so that . The viscous shallow water equations are then given as:


where is the acceleration of gravity, is the identity matrix, and is the dynamic dissipation coefficient to be defined shortly. In (1a), the coefficient defines whether viscosity is turned on () or off () in the continuity equation. We will return to this topic shortly.

3 Space and time discretization

High-order spectral element (SEM or CG, for continuous Galerkin) and discontinuous Galerkin (DG) approximations on quadrilateral elements are used to discretize Eq. (1).

The numerical model used in this paper is in fact the two-dimensional version of the NUMA model described in Giraldo et al. giraldoHesthavenWarburton2002 () and in Abdi and Giraldo abdiGiraldo2016 (). Furthermore, the model is derived from the model described in Kopera and Giraldo koperaGiraldo2013a (); koperaGiraldo2013b () which is here used for the shallow water equations. The solution is advanced in time using a fully implicit Runge-Kutta scheme (see §3.2).

3.1 Spectral element and discontinuous Galerkin approximations

We leave the details of the discretization to the work of Giraldo et al. giraldoHesthavenWarburton2002 () and Kopera and Giraldo koperaGiraldo2013a (); koperaGiraldo2013b (). Nonetheless, we introduce some notation that we are going to use later in the paper. To solve the shallow water equations by element-based Galerkin methods on a domain , we proceed by defining the weak form of (1) that we first recast in compact notation as


where is the transposed array of the solution variables and and are the flux and source terms.

In the case of spectral elements, the space discretization yields the semi-discrete matrix problem


where, for the global mass and differentiation matrices, and , we have that . The global matrices are obtained from their local counterparts, and , by direct stiffness summation, which maps the local degrees of freedom of an element to the corresponding global degrees of freedom in , and adds the element values in the global system. By construction, is diagonal (assuming inexact integration), with an obvious advantage if explicit time integration is used.

In the discontinuous Galerkin approximation, the problem at hand is solved only locally, and unlike the case of spectral elements, the flux integral that stems from the integration by-parts must be discretized as well. The element-wise counterpart of the matrix problem (3) is then written as:


where we obtain from the element boundary matrix, , and the element mass matrix, . Out of various possible choices for the definition of the numerical flux in Eq. (4), we adopted the Rusanov flux. The Laplace operator of viscosity is approximated using the Symmetric Interior Penalty method (SIP; the reader is referred to Arnold arnold1982 () for details on its definition).

3.2 Time integration

Equation (3) is integrated in time by an implicit Runge-Kutta (RK) scheme that corresponds to the implicit part of the implicit-explicit scheme used in giraldoEtAl2013 () (see also Butcher_1999 ()). The method coefficients in standard () tableaux form are the following


Scheme (5) is a three stage second order explicit-first-stage singly diagonally implicit RK (ESDIRK) scheme. This scheme has desirable accuracy and stability properties: all stages are second order accurate, it is stiffly accurate and -stable, and it is strong-stability-preserving (SSP) gottlieb2001strong () with SSP coefficient of 2. These properties allow us to take large time-steps with high accuracy as well as alleviate the instability issues associated with sharp solution gradients gottlieb2001strong (). The two-dimensional tests presented later in this paper showed to be the most demanding in terms of stability constraints. Method (5) allows us to gain up to one order of magnitude in terms of maximum admissible advective CFL when compared to an explicit method (explicit part of ARK3, Kennedy_2001 ()). In particular, the explicit four-stage Runge-Kutta solution of the solitary wave against one isolated obstacle described in §6.5 preserved stability for up to CFL=0.21 using both CG and DG approximations. Although we were not able to use arbitrarily large time-steps with the ESDIRK with the current implementation (we will address this issue in a future work), we resolved the same problem at CFL = 1.8. Schemes with a subset of these properties are employed by Kärnä et al. karnaEtAl2011 () and shown to be robust in this context. Method (5) used in this study satisfies all properties (-).

Computationally, at each of the two implicit stages we have to solve a nonlinear equation , where are the stage values, . We do so by using Newton iterations with a stopping criterion based on the relative decrease in the residual; that is, stop at iteration if . At each Newton iteration we have to solve a linear system , where is the Jacobian matrix of . We approximate the Jacobian using directional finite differences and iterate with the generalized minimal residual (GMRES) method, which is effectively a Jacobian-free Newton–Krylov method knoll2004jacobian (). The GMRES stopping criterion is also based on the relative residual. The first stage is explicit and equal to the last stage of the previous step, effectively making it a two-stage method, which saves some computational time.

4 The Dyn-SGS model for the shallow water equations

There are different ways to derive the viscous model described by Eq. (1) from the inviscid Saint-Venant equations. Analogous to our previous work on the large eddy simulation of stratified atmospheric flows marrasEtAl2015b (), the current model builds on the separation between grid resolved (indicated as for any quantity ) and unresolved (sub-grid) scales. The unresolved scales are modeled via the eddy viscosity scheme described in this paper. Given an element of order and with side lengths of comparable orders of magnitude, we define the following characteristic length (and hence filter width sagautBook ()):

The value of sets the size of the smallest resolvable scales in the same way as cut-off filters do in large eddy simulation models.

The application of this filter to the continuity equation (1a) results in the presence of an additional sub-grid term on the right-hand side. It is often debated whether diffusion should be applied to the continuity equation pasquettiGuermondPopoICOSAHOM2014 (); gerbeauEtAl2001 (); Guermond_Popov_2014 (); clearly, should the discrete viscous operator not be conservative, mass dissipation should not be used. However, by relying on spectral elements with integration by parts of the second-order diffusion operator, the discrete viscous operator is conservative, as shown in gubaEtAl2014 (). To get a sense of how necessary a stabilized continuity equation may be, we will show a few results for both conditions in §6.5.

Scale separation in the momentum equation yields a new equation that includes the gradient of the quantity


where the indicates sub-grid velocity. The coefficient is defined element-wise and is given as:






In (8, 9), indicates the spatially averaged value of the quantity at hand over the global domain , the norms at the denominator are used to preserve the physical dimension of the resulting equation, and are the residuals of the inviscid governing equations. At each time-step, the residuals are known. The presence of makes the artificial diffusion mathematically consistent111In the context of stabilization, mathematical consistency is the property of artificial diffusion such that it vanishes when the residual is zero.. The quantity in (8, 9) is the maximum wave speed.

Remark 1: physical dimensions

We would like to emphasize the necessity for the physically correct dimensions of diffusion. This is an important issue that is often underestimated and not accounted for in the design of artificial diffusion methods for stabilization purposes.

5 Wetting and drying interface

The dry regions are handled by adding an infinitesimal layer of water on the dry surfaces. We utilized a threshold value of 1e-3 m222One millimeter of water is physically negligible. for all the tests that we ran. Furthermore, the limiter of Xing et al. xingZhangShu2010 () is applied on the velocity and water depth. Even at high order, this approach is suitable for wetting and drying problems solved via both continuous and discontinuous Galerkin methods.

Albeit simple, the Xing-Zhang-Shu limiter works well. However, it is guaranteed to work if and only if the mean water height is positive. Xing et al. xingZhangShu2010 () recommend using the total variation bounding (TVB) limiter of Shu shu1987 (); cockburnShu1989 () to make the solution positive definite. We rely on artificial dissipation to achieve this as the TVB limiter requires a bound on the derivative of the solution, which we do not know a priori.

6 Numerical tests

We verify the correctness of our numerical model through five standard benchmarks in both one- and two-dimensions. In 1D, given the N-wave by Carrier et al. carrierEtAl2003 () that mimics a wave generated by an offshore submarine landslide, we compute the solution of the tsunami approaching a sloping beach. See §6.1. The numerical solution is evaluated against a set of tabulated surface height and momentum data provided by tsunamiWorkshopData (). The next 1D problem consists of the oscillation of a flat lake in a parabolic bowl gallardoEtAl2007 (); xingZhangShu2010 (). The existence of an analytic solution to this problem allows us to measure the accuracy of the numerical solutions. See §6.2. In addition, various 1D wetting and drying test cases were analyzed in the thesis of Hood karolineTHESIS2016 () to assess our model.

In 2D, we analyzed the results for the following tests. The vertical oscillation of a column of water in an axisymmetric paraboloid with analytic solution thacker1981 (). We describe this problem in §6.3 to test the inundation against the two-dimensional effects of a varying topography. The second two-dimensional test involves the flooding in a closed channel as presented by Gallardo et al. gallardoEtAl2007 () and by Xing and Zhang xingZhang2013 (). No exact solution is presented for this test, but the results can be compared against previous studies. This is shown in §6.4. We conclude the set of 2D problems with the simulation of the flow in a two-dimensional channel with an isolated obstacle synolakis1987 (). We use this test to gain insight into the practical – other than theoretical as stated in gerbeauEtAl2001 () – need to include viscosity by analyzing the energetics of the solutions. The test is run during a sufficiently long time to assess the robustness of the inundation and viscosity schemes to handle fast motion with interacting waves that impinge against a steep obstacle. See §6.5. The description and analysis of more tests can be found in the collection recently compiled by Delestre et al. delestreEtAl2014 ().

6.1 1D tsunami runup over a sloping beach

The runup of a long wave on a uniformly sloping beach was originally proposed at the third international workshop on long-wave runup models tsunamiWorkshopData (). The one-dimensional computational domain is defined as m. The dry initial beach is 500 m long. The initial waveform was defined by Carrier et al. in carrierEtAl2003 () for an m domain as:


with constants . To scale the wave to the m long domain used for the current test, the scaling factor is introduced and Eq. (10) is re-expressed with respect to and , with larger amplitudes . The initial wave is plotted in Fig. 1.

Figure 1: 1D tsunami runup over a sloping beach: initial N-wave.

The solutions at times s are plotted in Fig. 2 and are compared against the tabulated data available in tsunamiWorkshopData (). Fig. 2 shows that the effect of diffusion on the water surface solution is clearly negligible. This can be explained by looking at the structure and values of in Fig. 3. With a water surface that is smooth almost everywhere, the dynamic diffusion coefficient is so small that its effect becomes minimal. We will see later that this will not be the case in problems with a greater degree of irregularity of the surface.

We plot the space-time evolution of the shoreline in Fig. 4 where the wave elevation and total depth are plotted in the proximity of the shore. The dashed red curve represents the tabulated shoreline tsunamiWorkshopData (). By direct comparison with Carrier’s results carrierEtAl2003 (), the patterns of the water surface elevation () and total water depth are in full agreement.

Figure 2: 1D tsunami runup over a sloping beach. Computed and exact solutions at s. Left column: m (1250 elements of order ). Right column: m (2500 elements of order ). The computed inviscid and viscous (SGS) solutions appear almost perfectly superimposed. The problem is smooth almost everywhere so that the intensity of the dynamic dissipation is minimal. This translates into an almost imperceptible effect of diffusion on the solution. The structure of the dissipation coefficient is plotted in Fig. 3.
Figure 3: 1D tsunami runup over a sloping beach. Dynamic (red, dashed line) and water surface (blue, solid line) in the full domain. The effect of diffusion on the solution of Fig. 2 is minimal as the value of the coefficient is indeed very small. The solution is smooth almost everywhere, which is the reason for the very small values of the dynamic diffusion coefficient.
Figure 4: 1D tsunami runup over a sloping beach. Left: wave trajectory in the full 50 km long domain. Right: variation of the water depth in the proximity of the coast. The shoreline is at the interface between the white area (dry shore) and the color shading (water surface). The dashed red curve is the exact shoreline.

6.2 Grid convergence rate

To measure the convergence rate of the model, we compare the computed solutions against the analytic solution of the flow in a one-dimensional parabolic bowl gallardoEtAl2007 (); thacker1981 (). The parabolic topography is defined as:

where m and in m. The initial velocity is and the water surface begins to oscillate due to gravity only. The solution is computed using 8, 16, 32, 64, and 128 elements of order 4. The CG solution using 128 elements is plotted in Fig. 5, where we compare the inviscid solution against the viscous computation at s. The solution preserves its smoothness at all times as long as the resolution is sufficiently fine. At high resolution, it is evident that the CG approximation to this problem does not require stabilization. As the resolution is decreased by a factor of 4 (Fig. 6), the inviscid solution begins to develop oscillations although stability is not compromised. Unlike the inviscid one, the stabilized solution preserves the surface flatness although it moves slower due to a viscous source term in the momentum equation. This fact is also reflected in the computation of the normalized error norms plotted in Fig. 7. The analytic solution against which the error is calculated is defined for an inviscid flow. A difference between the viscous and inviscid computation is hence expected. The same observation applies to the CG and DG computations alike. The DG solutions to the same problem are plotted in Fig. 8 and 9.

Figure 5: 1D flow in a parabolic bowl. CG solution. Left column: inviscid. Right column: viscous. Computed water level (blue, solid line), exact solution (blue, dashed line with open squares), and (red, solid line) at different times using 128 elements of order 4. For visualization, is scaled by a factor of 1000. The solution is sufficiently smooth that the dynamic viscosity is negligible. On the other hand, the limiterwetting/drying scheme is responsible for the positivity of the water surface in the proximity of the moving boundary.
Figure 6: 1D flow in a parabolic bowl. Low resolution CG solutions using 16 elements of order 4. Left column: inviscid. Right column: viscous. Computed water level (blue, solid line), exact solution (blue, dashed line with open squares), and (red, solid line) at s. For visualization, is scaled by a factor of 1000.
Figure 7: Normalized error of water surface at s for CG (left) and DG (right). The and curves indicate the reference rates.
Figure 8: 1D flow in a parabolic bowl. Like Fig. 5 but using DG.
Figure 9: Like Fig. 6, but using DG.

6.3 2D oscillation in a symmetrical paraboloid

This test was defined by Thacker thacker1981 () who computed the exact solution of the problem. The paraboloid of revolution is defined as :

in the domain , with m. The initial condition is a reversed paraboloid (See top row of Fig. 11). The variation of the bed slope and the radial symmetry are a perfect test to assess the inundation scheme with two dimensional effects. The inviscid CG and DG solutions are plotted in Figs. 10 and 11, respectively, between and s. The radial symmetry of the solution is preserved throughout the simulation, and no spurious modes can be observed in the proximity of the moving shoreline. Although it is not fully visible from the plots, we observed a small discrepancy between the exact and computed solutions that occurs at the very center of the bowl at the latest time. More specifically, the numerical solution tends to move slightly slower than the exact solution along the bowl centerline. This occurs for CG and DG alike.

Figure 10: 2D oscillation in a paraboloid. Inviscid CG solution. Left: view. Right: view. From top to bottom: . Solution with elements of order 4.
Figure 11: 2D oscillation in a paraboloid. Inviscid DG solution. Left: view. Right: view. From top to bottom: . Solution with elements of order 4.

6.4 2D flooding problem in a channel with three mounds

This test kawaharaUmetsu1986 (); brufauEtAl2002 () was used by Xing and Zhang xingZhang2013 () and Gallardo et al. gallardoEtAl2007 () to assess, on unstructured grids, the same positivity-preserving method that we are applying here xingZhangShu2010 (). We reproduce their results here to verify its correct implementation in our code. The three mounds in the channel , are defined by the function where


The flooding is triggered by a dam break at s. The evolution of the flooding from to s is shown in Fig. 12. We computed the solution using a relatively coarse grid made of elements of order 4. The results are in agreement with gallardoEtAl2007 () and xingZhang2013 (): the flow separation is well resolved and occurs after 10 seconds from the dam breaking. After 15 seconds, the first water front has reached the back wall and is fully reflecting back by 20 seconds. By 40 seconds, the flow has almost reached full rest.

Figure 12: Dam break problem in a closed channel with three mounds. Viscous CG solution at s. Solution computed using elements of order 4 for the domain .

6.5 2D solitary wave runup and run-down on a circular island

A solitary wave runup on a circular island was studied in briggsEtAl1995 () at the Waterways Experiment Station of the US Army Corps of Engineers. In this example, the initial wave is modeled via the following analytic definition by Synolakis synolakis1987 ():


where m is the wave amplitude, m, m is the initial still water level, and

The island is a cone given as

where , m, and is centered at m. The cone is installed on a flat bathymetry. The fluid is confined within four solid walls.

To understand how the proposed diffusion and numerical approximation depend on the grid, we ran the simulation at the four resolutions m. Fig. 13 shows that the stabilized DG solution is converging to the same solution and the main features of the interacting waves are reproduced almost equally across the four grids. Certainly, the 40 cm grid spacing is the most dissipative, although it is encouraging to see how the important features resolved at 5 cm are still well represented on the coarsest grid. The same observation applies to the CG solution (plot not shown).

Figure 13: DG solution of water depth for the single-hill configuration. Results obtained with the four grid resolutions m (indicated in the figures) using -order elements. The color bar is cut at 0.25 to preserve the visibility of the smallest features. The dark blue coloring within the region of the cone corresponds to the water depth equal to the threshold water layer of 1e-3 m.

In Figs. 14-17 we plot the projection of the 2D solution on the plane . The spurious modes that characterize the water surface in the proximity of the sharpest wave front are fully removed by Dyn-SGS (Fig. 14) without weakening the front sharpness. We observe the same in the velocity field plotted in Fig. 15. This is in full agreement with the application of Dyn-SGS to non-linear wave problems with strong discontinuities, as previously shown in (marrasEtAl2015b, , Figs. 16, 17) for the solution of the Burgers’ equation. We briefly mentioned above how DG already has built-in dissipation. This is clearly visible from the plot of Fig. 16; the unstabilized DG solution shows no oscillations except for, at most, some minimal under- and over-shooting. This implies that its residual is so small that the effect of SGS eddy viscosity reduces to a minimum value. This explains why the inviscid and viscous DG solutions look similar.

It is important to notice the same space distribution of the CG and DG velocity waves shown in Figs. 15 and 17. This is a visual confirmation that the two solutions present the same dispersion properties, which is to be expected from two analogous numerical approximations to the same problem. This is indicative of a correct implementation of the unified CG/DG framework.

Figure 14: view of the 2D CG solution of water depth for the single-hill configuration. Inviscid (left) against viscous solution using SGS (right). Solutions obtained using -order elements. The dark blue coloring within the region of the cone corresponds to the water depth equal to the threshold water layer of 1e-3 m.
Figure 15: Like Fig. 14, but for -velocity component. Inviscid (left) against viscous solution using SGS (right).
Figure 16: view of the 2D DG solution of water depth for the single-hill configuration. Inviscid (left) against viscous solution using SGS (right). Solution obtained using -order elements. The dark blue coloring within the region of the cone corresponds to the water depth equal to the threshold water layer of 1e-3 m. This plot shows the power of DG. Without SGS it still almost captures the bore sharply.
Figure 17: Like Fig. 16, but for -velocity component. Inviscid (left) against viscous solution using SGS (right). A comparison of this figure with Fig. 15 shows the same dispersion properties of the CG and DG solutions.

The solution of the top plot in Fig. 18, was computed with stabilization applied to the continuity ( in Eq. (1)) and momentum equations. When compared against the un-stabilized solution (bottom plot in the same figure), we notice that the features of the fronts of the interacting waves are fully preserved. Furthermore, the fronts are not excessively smeared out as the high frequency modes are removed. The effect of momentum stabilization on the velocity field is shown in Fig. 19.

Figure 18: CG solution for the single-hill configuration. Instantaneous perspective view of the unstabilized (top) and stabilized (bottom) water surface for cm using -order elements. The high frequency instabilities are removed by Dyn-SGS without compromising the overall sharpness of the interacting waves. Both solutions are characterized by the same wave features at all scales.
Figure 19: Like Fig. 18, but surface of horizontal velocity, . For best view of the velocity surface, the view angle is different from the one of Fig. 18.

The instantaneous energy spectra of the viscous and inviscid CG and DG solutions at s are plotted in Fig. 20. The difference between the CG and DG curves is striking. The viscous and inviscid DG spectra overlap almost fully and show approximately the same decay across the entire spectrum, from a -5/3 slope in the inertial sub-range to a -3 slope in the dissipation wave numbers (refer to boffettaEcke2012 () for a review on two-dimensional flows and their energetics). This is only true as long as the resolution is not too coarse, especially so in the case of CG. At very coarse resolutions ( m), neither CG or DG can avoid energy from building up in the highest modes unless artificial diffusion is used. The inherent dissipation of DG is no longer sufficient.

Figure 20: Instantaneous 1D energy spectra of the single hill problem of Fig. 13 at s. Left: CG with and without viscosity. Right: DG with and without viscosity. From top to bottom the resolution decreases.

We stated above that is only active where the equation residuals (i.e. gradients) are important. In the case of water waves, this occurs in the proximity of the wave fronts. In Fig. 21, we plot to show its spatial structure and its evolution between and seconds. This plots clearly show how diffusion is equally zero away from the fronts and only activates where really necessary. It may not be so obvious to achieve this by using an artificial diffusion that is not residual-based. To provide a visual correlation between and the wave features, in Fig. 22 we plot the stabilized spectral element solution of the water surface at a grid resolution m.

Figure 21: Time evolution of from s to s. The colorscale ranges between 0 (blue) and 5 (red) ms. The plotted domain is m.
Figure 22: Time evolution of from s to s. The color scale ranges between 0.25 (blue) and 0.5 (red) m. The plotted domain is m.

7 Conclusions

We presented the numerical solution of the shallow water equations via continuous and discontinuous Galerkin (CG/DG) methods to model problems involving inundation. Careful handling of the transition between dry and wet surfaces is particularly challenging for high-order methods. The most adopted solution in these cases consists of lowering the approximation order in either the transition zone alone, or in the whole domain. By extending to CG a simple slope limiter originally designed for DG xingZhangShu2010 (), and by combining it with a thin water layer always present in dry regions, we showed that the nominal high-order accuracy of the underlying space approximations was preserved and that the solution smoothness in the proximity of moving boundaries was maintained in one- and two-dimensions. The simplicity of this approach is effective for CG as well as DG. More important, we demonstrated its effectiveness using higher order elements. This limiter, unless embedded with a positivity-preserving scheme for water height, does not prevent the high-order solution from triggering unphysical Gibbs oscillations in the proximity of strong gradients. To overcome this problem, we presented a dynamically adaptive dissipation based on a residual-based sub-grid scale eddy viscosity model (Dyn-SGS). By numerical examples, we demonstrated the following properties of this model when solving the shallow water equations: It removes the Gibbs oscillations that form in the proximity of sharp wave fronts while preserving the overall accuracy and sharpness of the solution everywhere else. This is possible thanks to the residual-based definition of the dynamic diffusion coefficient. For coarse grids, it prevents energy from building up at small wave-numbers; this is very important to preserve numerical stability in the flow regimes we are interested in. The model has no user tunable parameter, which is of great advantage when the model is to be used by an external user. It is important to underline that the natural, built-in dissipation of DG be may be large enough that the contribution of Dyn-SGS becomes irrelevant. When this happens, the dynamic dissipation detects it from the residual, and hence limits its own strength. Finally, a three stage, second order explicit-first-stage, singly diagonally implicit Runge-Kutta (ESDIRK) time integration scheme was implemented to overcome the small time-step restriction incurred by high-order Galerkin approximations.

8 Acknowledgments

The authors would like to acknowledge the contribution of Haley Lane, who implemented the one-dimensional version of the wetting and drying algorithm used in this work. The authors would also like to acknowledge Karoline Hood who tested the correctness of the implicit solver in her NPS Master’s thesis (karolineTHESIS2016 ()). The authors are also thankful to Oliver Fringer and the Rogers for discussions regarding coastal flows, and to Stephen R. Guimond for providing his MATLAB functions to compute the energy spectra. FXG acknowledges the support of the ONR Computational Mathematics program, and FXG and EMC acknowledge the support of AFOSR Computational Mathematics.


  • (1) The third international workshop on long-wave runup models: Benchmark problem #1: Tsunami runup onto a plane beach. http://isec.nacse.org/workshop/2004_cornell/bmark1.html. 2004.
  • (2) D. Abdi, L. C. Wilcox, T. Warburton, and F. X. Giraldo. A GPU accelerated continuous and discontinuous Galerkin non-hydrostatic atmospheric model. Submitted. Pre-print available at www.researchgate.net/publication/
    , 2016.
  • (3) D. S. Abdi and F. X. Giraldo. Efficient construction of unified continuous and discontinuous Galerkin formulations for the 3D Euler equations. J. Comput. Phys., 320:46–68, 2016.
  • (4) M. Ainsworth, P. Monk, and W. Muniz. Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation. J. Sci. Comput., 27:5–40, 2006.
  • (5) D. Arnold. An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal., 19:742–760, 1982.
  • (6) N. Beisiegel and B. Behrens. Quasi-nodal third-order Bernstein polynomioals in a discontinuous Galerkin model for flooding and drying. Environ. Earth Sci., 74:7275–7284, 2015.
  • (7) G. Boffetta and E. Ecke. Two-dimensional turbulence. Annu. Rev. Fluid Mech., 44:427–451, 2012.
  • (8) J. P. Boyd. The erfc-log filter and the asymptotics of the Euler and Vandeven sequence accelerations. In A.V. Ilin, L.R. Scott (Eds.), Proceedings of the Third International Conference on Spectral and High Order Methods, Houston Journal of Mathematics, pages 267–276, 1996.
  • (9) M. J. Briggs, C. E. Synolakis, G. S. Harkins, and D. R. Green. Laboratory experiments of tsunami runup on a circular islad. Pageoph., 144:569–593, 1995.
  • (10) P. Brufau, M. E. Vázquez-Cendón, and P. García-Navarro. A numerical model for the flooding and drying of irregular domains. Int. J. Numer. Methods Fluids, 39:247–275, 2002.
  • (11) S. Bunya, E. J. Kubatko, J. J. Westerink, and C. Dawson. A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations. Comput. Methods Appl. Mech. Engrg., 198:1548–1562, 2009.
  • (12) J.C. Butcher and P. Chartier. The effective order of singly-implicit Runge-Kutta methods. Numerical Algorithms, 20(4):269–284, 1999.
  • (13) G. F. Carrier, T. T. Wu, and H. Yeh. Tsunami run-up and draw-down on a plane beach. J. Fluid Mech., 475:79–99, 2003.
  • (14) B. Cockburn and C W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Math. Comput., 52:411–435, 1989.
  • (15) C Dawson and V. Aizinger. A discontinuous Galerkin method for three-dimensional shallow water equations. J. Sci. Comp., 22:245–267, 2995.
  • (16) A. J. C. de Saint-Venant. Théorie du mouvement non-permanent de eaux, avec application aux crues de rivières et à l’introduction des marées dans leur lit. C. R. Acad. Sc. Paris, 73:147–154, 1871.
  • (17) O. Delestre, C. Lucas, P. A. Ksinant, F. Darboux, C. Laguerre, T. N. T. Vo, F. James, and S. Corier. SWASHES: a compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies. arXiv, 1110.0288v6 [math.NA], 2014.
  • (18) D. Drikakis, H. Marco, F. F. Grinstein, C. R. DeVore, C. Fureby, M. Liefvendahl, and D. L. Youngs. Numerics for ILES: limiting algorithms. In F. F. Grinstein, L. G. Margolin, and Rider W. J., editors, Implicit Large-Eddy Simulation: Computing Turbulent Flow Dynamics, pages 94–129. Cambridge University Press, 2007.
  • (19) J. M. Gallardo, C. Parés, and M. Castro. On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas. J. Comput. Phys., 227:574–601, 2007.
  • (20) R. Gandham, D. Medina, and T. Warburton. GPU accelerated discontinuous galerkin methods for shallow water equations. Commun. Comput. Phys., 18:37–64, 2015.
  • (21) J. F. Gerbeau and B. Perthame. Derivation of the viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete Contin. Dyn. Syst. Ser. B, 1:89–102, 2001.
  • (22) F. X. Giraldo. A spectral element shallow water model on spherical geodesic grids. Int. J. Num. Meth. Fluids, 35:869–901, 2001.
  • (23) F. X. Giraldo, J. S. Hesthaven, and T. Warburton. Nodal high-order discontinuous Galerkin methods for spherical shallow water equations. J. Comput. Phys., 181:499–525, 2002.
  • (24) F X. Giraldo, J F.. Kelly, and E. Constantinescu. Implicit-explicit formulations of a three-dimensional Nonhydrostatic Unified Model of the Atmosphere (NUMA). SIAM J. Sci. Comput., 35:1162–1194, 2013.
  • (25) F. X. Giraldo and M. Restelli. High-order semi-implicit time-integrators for a triangular discontinuous Galerkin oceanic shallow water model. Int. J. Numer. Methods Fluids, 63:1077–1102, 2010.
  • (26) S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. SIAM review, 43(1):89–112, 2001.
  • (27) R. Gourgue, O. Comblen, J. Lambrechts, T. Kärnä, V. Legat, and E. Deleersnijder. A flux-limiting wetting-drying method for finite-element shallow-water models, with application to the Scheldt Estuary. Adv. Water Resourc., 32:1726–1739, 2009.
  • (28) O. Guba, M. A. Taylor, P. A. Ullrich, J. R. Overfelt, and M. N. Levy. The spectral element method on variable resolution grids: evaluating grid sensitivity and resolution-aware numerical viscosity. Geosci. Model Dev. Discuss., 7:4081–4117, 2014.
  • (29) J L. Guermond and B. Popov. Viscous regularization of the Euler equations and entropy principles. SIAM J. Appl. Math., 74(2):284–305, 2014.
  • (30) E. Hendricks, M. Kopera, and F. X. Giraldo. Evaluation of the utility of static and adaptive mesh refinement for idealized tropical cyclone problems in a spectral element shallow water model. Mon. Wea. Rev. (in press), 2016.
  • (31) K. Hood. Modeling storm surge using discontinuous Galerkin methods. Masters Thesis, Dept. of Applied Maths, Naval Postgraduate School, Monterey, CA, U.S.A., 2016.
  • (32) M. Iskandarani, D .B. Haidvogel, and J. P. Boyd. A staggered spectral element models with application to the oceanic shallow water equations. Int. J. Numer. Methods Fluids, 20:393–414, 1995.
  • (33) T. Kärnä, B. de Brye, O. Gourgue, J. Lambrechts, R. Comblen, V. Legat, and E. Deleersnijder. A fully implicit wetting-drying method for DG-FEM shallow water models, with an application to the Scheldt Estuary. Comp. Methods Appl. Mech. Engrg., 200:509–524, 2011.
  • (34) M. Kawahara and T. Umetsu. Finite element method for moving boundary problems in river flow. Int. J. Numer. Meth. Fluids, 6:365–386, 1986.
  • (35) C.A. Kennedy and M.H. Carpenter. Additive Runge-Kutta schemes for convection-diffusion-reaction equations. Appl. Numer. Math., 44(1-2):139–181, 2003.
  • (36) D.A. Knoll and D.E. Keyes. Jacobian-free Newton–Krylov methods: a survey of approaches and applications. Journal of Computational Physics, 193(2):357–397, 2004.
  • (37) M. A. Kopera and F. X. Giraldo. Analysis of adaptive mesh refinement for IMEX discontinuous Galerkin solutions of the compressible Euler equations with application to atmospheric simulations. J. Comput. Phys., 92-117:275, 2014.
  • (38) M. A. Kopera and F. X. Giraldo. Mass conservation of the unified continuous and discontinuous element-based galerkin methods on dynamically adaptive grids with application to atmospheric simulations. J. Comput. Phys., 297:90–103, 2015.
  • (39) E. J. Kubatko, J. J. Westerink, and C. Dawson. hp discontinuous Galerkin methods for advection dominated problems in shallow water flow. Comp. Methods Appl. Mech. Engrg., 196:437–451, 2006.
  • (40) D. K. Lilly. On the numerical simulation of buoyant convection. Tellus, 14:148–172, 1962.
  • (41) H. Ma. A spectral element basin model for the shallow water equations. J. Comput. Phys., 109:133–149, 1993.
  • (42) S. Marras, M. Nazarov, and F. X. Giraldo. Stabilized high-order Galerkin methods based on a parameter-free dynamic SGS model for LES. J. Comput. Phys., 301:77–101, 2015.
  • (43) C. Michoski, C. Dawson, E. J. Kubatko, D. Wirasaet, S. Brus, and J. J. Westerink. A comparison of artificial viscosity, limiters, and filters, for high order discontinuous Galerkin solutions in nonlinear settings. J. Sci. Comput., 66:406–434, 2016.
  • (44) A. Müller, M. A. Kopera, S. Marras, L .C. Wilcox, T. Isaac, and F. X. Giraldo. Strong Scaling for Numerical Weather Prediction at Petascale with the Atmospheric Model NUMA. Submitted to the IEEE International Parallel and Distributed Processing Symposium. Pre-print available at www.researchgate.net/publication282887372_
    , 2016.
  • (45) R D. Nair, S J. Thomas, and R D. Loft. A discontinuous Galerkin global shallow water model. Mon. Wea. Rev., 133, 2007.
  • (46) M. Nazarov and J. Hoffman. Residual-based artificial viscosity for simulation of turbulent compressible flow using adaptive finite element methods. Int. J. Numer. Methods Fluids, 71:339–357, 2013.
  • (47) R. Pasquetti, J.L. Guermond, and B. Popov. Stabilized spectral element approximation of the Saint Venant system using the entropy viscosity technique. In Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014, Lecture Notes in Computational Science and Engineering, pages 397–404, 2015.
  • (48) C. Phan Van, E. Deleersnijder, D. Bousmar, and S. Frazão Soares. Simulation of flow in compound open-channel using a discontinuous Galerkin finite-element method with Smagorinsky turbulence closure. J. Hydro-env. Research, 8:396–409, 2014.
  • (49) N. Rakowsky, A. Androsov, A. Fuchs, S. Harig, A. Immerz, S. Danilov, W. Hiller, and J. Schröter. Operational tsunami modelling with TsunAWI - recent developments and applications. Nat. Haz. Earth Sys. Sci., 13:1629–1642, 2013.
  • (50) W. J. Rider. The relationship of MPDATA to other high-resolution methods. Int. J. Numer. Methods Fluids, 50:1145–1158, 2006.
  • (51) P. Sagaut. Large Eddy Simulation for Incompressible flows. Springer, 3rd edition, 2006.
  • (52) C W. Shu. TVB uniformly high-order schemes for conservation laws. Math. Comput., 49:105–121, 1987.
  • (53) J. Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiement. Mon. Wea. Rev., 91:99–164, 1963.
  • (54) P. K. Smolarkiewicz. A simple positive definite advection scheme with small implicit diffusion. Mon. Wea. Rev., 111:479–486, 1983.
  • (55) P. K. Smolarkiewicz. A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. J. Comput. Phys., 54:325–362, 1984.
  • (56) C. E. Synolakis. The runup of solitary waves. J. Fluid Mech., 185:523–545, 1987.
  • (57) M. Taylor, J. Tribbia, and M. Iskandarani. The spectral element method for the shallow water equations on the sphere. J. Comput. Phys., 130:92–108, 1997.
  • (58) W. C. Thacker. Some exact solutions to the nonlinear shallow-water wave equations. J. Fluid Mech., 107:499–508, 1981.
  • (59) A. Uranga, P.-O. Persson, M. Drela, and J. Peraire. Implicit Large Eddy Simulation of transition to turbulence at low Reynolds numbers using a discontinuous Galerkin method. Int. J. Numer. Methods Engrg., 87:232–261, 2011.
  • (60) H. Vandeven. Family of spectral filters for discontinuous problems. J. Sci. Comp., 6:159–192, 1991.
  • (61) S. Vater, N. Beisiegel, and J. Behrens. A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows with wetting and drying: one dimensional case. Adv. Water Res., 85:1–13, 2015.
  • (62) Y. Xing and X. Zhang. Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes. J. Sci. Comput., 57:19–41, 2013.
  • (63) Y. Xing, X. Zhang, and C. W. Shu. Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations. Adv. Water Res., 33:1476–1493, 2010.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description