On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows

We present a new unconditionally positivity‐preserving (PP) implicit time integration method for the DG scheme applied to shallow water flows. This novel time discretization enhances the currently used PP DG schemes, because in the majority of previous work, explicit time stepping is implemented to deal with wetting and drying. However, for explicit time integration, linear stability requires very small time steps. Especially for locally refined grids, the stiff system resulting from space discretization makes implicit or partially implicit time stepping absolutely necessary. As implicit schemes require a lot of computational time solving large systems of nonlinear equations, a much larger time step is necessary to beat explicit time stepping in terms of CPU time. Unfortunately, the current PP implicit schemes are subject to time step restrictions due to a so‐called strong stability preserving constraint. In this work, we hence give a novel approach to positivity preservation including its theoretical background. The new technique is based on the so‐called Patankar trick and guarantees non‐negativity of the water height for any time step size while still preserving conservativity. In the DG context, we prove consistency of the discretization as well as a truncation error of the third order away from the wet–dry transition. Because of the proposed modification, the implicit scheme can take full advantage of larger time steps and is able to beat explicit time stepping in terms of CPU time. The performance and accuracy of this new method are studied for several classical test cases. Copyright © 2014 John Wiley & Sons, Ltd.     

A dissipation‐free numerical method to solve one‐dimensional hyperbolic flow equations

In this paper, a numerical method to capture the shock wave propagation in 1‐dimensional fluid flow problems with 0 numerical dissipation is presented. Instead of using a traditional discrete grid, the new numerical method is built on a range‐discrete grid, which is obtained by a direct subdivision of values around the shock area. The range discrete grid consists of 2 types: continuous points and shock points. Numerical solution is achieved by tracking characteristics and shocks for the movements of continuous and shock points, respectively. Shocks can be generated or eliminated when triggering entropy conditions in a marking step. The method is conservative and total variation diminishing. We apply this new method to several examples, including solving Burgers equation for aerodynamics, Buckley‐Leverett equation for fractional flow in porous media, and the classical traffic flow. The solutions were verified against analytical solutions under simple conditions. Comparisons with several other traditional methods showed that the new method achieves a higher accuracy in capturing the shock while using much less grid number. The new method can serve as a fast tool to assess the shock wave propagation in various flow problems with good accuracy.     

A reduction method to quasilinear hyperbolic systems of multicomponent field PDEs with application to wave interaction

A reduction method is worked out for determining a class of exact solutions with inherent wave features to quasilinear hyperbolic homogeneous systems of N>2 first-order autonomous PDEs. A crucial point of the present approach is that in the process the original set of field equations induces the hyperbolicity of an auxiliary 2×2 subsystem and connection between the respective characteristic velocities can be established. The integration of this auxiliary subsystem via the hodograph method and through the use of the Riemann invariants provides the searched solutions to the full governing system. These solutions also represent invariant solutions associated with groups of translation of space/time coordinates and involving arbitrary functions that can be used for studying non-linear wave interaction. Within such a theoretical framework the two-dimensional motion of an adiabatic fluid is considered. For appropriate model pressure-entropy-density laws, we determine a solution to the governing system of equations which describes in the 2+1 space two non-linear waves which were initiated as plane waves, interact strongly on colliding but emerge with unaffected profile from the interaction region. These model material laws include the classical pressure-entropy-density law which is usually adopted for a polytropic fluid.     

An eigenvector‐based linear reconstruction scheme for the shallow‐water equations on two‐dimensional unstructured meshes

This paper presents a new approach to MUSCL reconstruction for solving the shallow‐water equations on two‐dimensional unstructured meshes. The approach takes advantage of the particular structure of the shallow‐water equations. Indeed, their hyperbolic nature allows the flow variables to be expressed as a linear combination of the eigenvectors of the system. The particularity of the shallow‐water equations is that the coefficients of this combination only depend upon the water depth. Reconstructing only the water depth with second‐order accuracy and using only a first‐order reconstruction for the flow velocity proves to be as accurate as the classical MUSCL approach. The method also appears to be more robust in cases with very strong depth gradients such as the propagation of a wave on a dry bed. Since only one reconstruction is needed (against three reconstructions in the MUSCL approach) the EVR method is shown to be 1.4–5 times as fast as the classical MUSCL scheme, depending on the computational application. Copyright © 2006 John Wiley & Sons, Ltd.     

Fifth‐order Hermitian schemes for computational linear aeroacoustics

We develop a class of fifth‐order methods to solve linear acoustics and/or aeroacoustics. Based on local Hermite polynomials, we investigate three competing strategies for solving hyperbolic linear problems with a fifth‐order accuracy. A one‐dimensional (1D) analysis in the Fourier series makes it possible to classify these possibilities. Then, numerical computations based on the 1D scalar advection equation support two possibilities in order to update the discrete variable and its first and second derivatives: the first one uses a procedure similar to that of Cauchy–Kovaleskaya (the ‘Δ‐P5 scheme’); the second one relies on a semi‐discrete form and evolves in time the discrete unknowns by using a five‐stage Runge–Kutta method (the ‘RGK‐P5 scheme’). Although the RGK‐P5 scheme shares the same local spatial interpolator with the Δ‐P5 scheme, it is algebraically simpler. However, it is shown numerically that its loss of compactness reduces its domain of stability. Both schemes are then extended to bi‐dimensional acoustics and aeroacoustics. Following the methodology validated in (J. Comput. Phys. 2005; 210 :133–170; J. Comput. Phys. 2006; 217 :530–562), we build an algorithm in three stages in order to optimize the procedure of discretization. In the ‘reconstruction stage’, we define a fifth‐order local spatial interpolator based on an upwind stencil. In the ‘decomposition stage’, we decompose the time derivatives into simple wave contributions. In the ‘evolution stage’, we use these fluctuations to update either by a Cauchy–Kovaleskaya procedure or by a five‐stage Runge–Kutta algorithm, the discrete variable and its derivatives. In this way, depending on the configuration of the ‘evolution stage’, two fifth‐order upwind Hermitian schemes are constructed. The effectiveness and the exactitude of both schemes are checked by their applications to several 2D problems in acoustics and aeroacoustics. In this aim, we compare the computational cost and the computation memory requirement for each solution. The RGK‐P5 appears as the best compromise between simplicity and accuracy, while the Δ‐P5 scheme is more accurate and less CPU time consuming, despite a greater algebraic complexity. Copyright © 2007 John Wiley & Sons, Ltd.     

A particle‐in‐cell method for the solution of two‐layer shallow‐water equations

A particle‐in‐cell (PIC) numerical method developed for the study of shallow‐water dynamics, when the moving fluid layer is laterally confined by the intersection of its top and bottom surfaces, is described. The effect of ambient rotation is included for application to geophysical fluids, particularly open‐ocean buoyant vortices in which the underlying density interface outcrops to the surface around the rim of the vortex. Extensions to include the dynamical effect of a second moving layer (baroclinicity) and the presence of a lateral rigid boundary (sidewall) are also described. Although the method was developed for oceanographic investigations, applications to other fluid mechanics problems would be straightforward. Copyright © 2000 John Wiley & Sons, Ltd.     

High‐order ADER‐WENO ALE schemes on unstructured triangular meshes—application of several node solvers to hydrodynamics and magnetohydrodynamics

In this paper, we present a class of high‐order accurate cell‐centered arbitrary Lagrangian–Eulerian (ALE) one‐step ADER weighted essentially non‐oscillatory (WENO) finite volume schemes for the solution of nonlinear hyperbolic conservation laws on two‐dimensional unstructured triangular meshes. High order of accuracy in space is achieved by a WENO reconstruction algorithm, while a local space–time Galerkin predictor allows the schemes to be high order accurate also in time by using an element‐local weak formulation of the governing PDE on moving meshes. The mesh motion can be computed by choosing among three different node solvers, which are for the first time compared with each other in this article: the node velocity may be obtained either (i) as an arithmetic average among the states surrounding the node, as suggested by Cheng and Shu, or (ii) as a solution of multiple one‐dimensional half‐Riemann problems around a vertex, as suggested by Maire, or (iii) by solving approximately a multidimensional Riemann problem around each vertex of the mesh using the genuinely multidimensional Harten–Lax–van Leer Riemann solver recently proposed by Balsara et al. Once the vertex velocity and thus the new node location have been determined by the node solver, the local mesh motion is then constructed by straight edges connecting the vertex positions at the old time level tⁿ with the new ones at the next time level t^{n + 1}. If necessary, a rezoning step can be introduced here to overcome mesh tangling or highly deformed elements. The final ALE finite volume scheme is based directly on a space–time conservation formulation of the governing PDE system, which therefore makes an additional remapping stage unnecessary, as the ALE fluxes already properly take into account the rezoned geometry. In this sense, our scheme falls into the category of direct ALE methods. Furthermore, the geometric conservation law is satisfied by the scheme by construction. We apply the high‐order algorithm presented in this paper to the Euler equations of compressible gas dynamics as well as to the ideal classical and relativistic magnetohydrodynamic equations. We show numerical convergence results up to fifth order of accuracy in space and time together with some classical numerical test problems for each hyperbolic system under consideration. Copyright © 2014 John Wiley & Sons, Ltd.     

Resolution to a three‐dimensional physical oceanographic problem using the non‐linear Galerkin method

This paper presents the numerical results concerning the adaptation of the non‐linear Galerkin method to three‐dimensional geophysical fluid equations. This method was developed by Marion and Temam to solve the Navier–Stokes two‐dimensional equations. It allows a substantial decrease in calculation costs due to the application of an appropriate treatment to each mode based on its position in the spectrum. The large scales involved in the study of geophysical flow require that the earth's rotational effects and the existence of a high degree of stratification be taken into account. These phenomena play an important role in the distribution of the energy spectrum. It is shown here that the non‐linear Galerkin method is very well‐suited to the treatment of these phenomena. First, the method for the particular situation of a rigid‐lid with a flat bottom is validated, for which the functional basis used is particularly well‐adapted. Then the more general case of a domain exhibiting variable bathymetry is presented, which necessitates the use of the transformation σ, thus providing a study domain with a cylindrical configuration. Copyright © 1999 John Wiley & Sons, Ltd.     

The numerical study of roll‐waves in inclined open channels and solitary wave run‐up

In this paper, we introduce a finite‐volume kinetic BGK scheme and its applications to the study of roll and solitary waves. The current scheme is based on the numerical solution of the gas‐kinetic Bhatnagar–Gross–Krook model in the flux evaluation across each cell interface. An intrinsic connection between the BGK model and time‐dependent, non‐linear, non‐homogeneous shallow‐water equations enables us to solve shallow‐water equations automatically with our kinetic scheme. The analytical solution, experimental measurements, and numerical calculations for problems associated with roll‐waves down an inclined open channel and solitary waves incident on a sloped beach are also presented. Copyright © 2005 John Wiley & Sons, Ltd.     

A particle finite element method applied to long wave run‐up

This paper presents a Lagrangian–Eulerian finite element formulation for solving fluid dynamics problems with moving boundaries and employs the method to long wave run‐up. The method is based on a set of Lagrangian particles which serve as moving nodes for the finite element mesh. Nodes at the moving shoreline are identified by the alpha shape concept which utilizes the distance from neighbouring nodes in different directions. An efficient triangulation technique is then used for the mesh generation at each time step. In order to validate the numerical method the code has been compared with analytical solutions and a preexisting finite difference model. The main focus of our investigation is to assess the numerical method through simulations of three‐dimensional dam break and long wave run‐up on curved beaches. Particularly the method is put to test for cases where different shoreline segments connect and produce a computational domain surrounding dry regions. Copyright © 2006 John Wiley & Sons, Ltd.     

Hybrid finite‐volume finite‐difference scheme for the solution of Boussinesq equations

A hybrid scheme composed of finite‐volume and finite‐difference methods is introduced for the solution of the Boussinesq equations. While the finite‐volume method with a Riemann solver is applied to the conservative part of the equations, the higher‐order Boussinesq terms are discretized using the finite‐difference scheme. Fourth‐order accuracy in space for the finite‐volume solution is achieved using the MUSCL‐TVD scheme. Within this, four limiters have been tested, of which van‐Leer limiter is found to be the most suitable. The Adams–Basforth third‐order predictor and Adams–Moulton fourth‐order corrector methods are used to obtain fourth‐order accuracy in time. A recently introduced surface gradient technique is employed for the treatment of the bottom slope. A new model ‘HYWAVE’, based on this hybrid solution, has been applied to a number of wave propagation examples, most of which are taken from previous studies. Examples include sinusoidal waves and bi‐chromatic wave propagation in deep water, sinusoidal wave propagation in shallow water and sinusoidal wave propagation from deep to shallow water demonstrating the linear shoaling properties of the model. Finally, sinusoidal wave propagation over a bar is simulated. The results are in good agreement with the theoretical expectations and published experimental results. Copyright © 2005 John Wiley & Sons, Ltd.     

Implicit schemes with large time step for non‐linear equations: application to river flow hydraulics

In this work, first‐order upwind implicit schemes are considered. The traditional tridiagonal scheme is rewritten as a sum of two bidiagonal schemes in order to produce a simpler method better suited for unsteady transcritical flows. On the other hand, the origin of the instabilities associated to the use of upwind implicit methods for shock propagations is identified and a new stability condition for non‐linear problems is proposed. This modification produces a robust, simple and accurate upwind semi‐explicit scheme suitable for discontinuous flows with high Courant–Friedrichs–Lewy (CFL) numbers. The discretization at the boundaries is based on the condition of global mass conservation thus enabling a fully conservative solution for all kind of boundary conditions. The performance of the proposed technique will be shown in the solution of the inviscid Burgers' equation, in an ideal dambreak test case, in some steady open channel flow test cases with analytical solution and in a realistic flood routing problem, where stable and accurate solutions will be presented using CFL values up to 100. Copyright © 2004 John Wiley & Sons, Ltd.     

A dual‐weighted trust‐region adaptive POD 4D‐VAR applied to a finite‐element shallow‐water equations model

We consider a limited‐area finite‐element discretization of the shallow‐water equations model. Our purpose in this paper is to solve an inverse problem for the above model controlling its initial conditions in presence of observations being assimilated in a time interval (window of assimilation). We then attempt to obtain a reduced‐order model of the above inverse problem, based on proper orthogonal decomposition (POD), referred to as POD 4‐D VAR. Different approaches of POD implementation of the reduced inverse problem are compared, including a dual‐weighed method for snapshot selection coupled with a trust‐region POD approach. Numerical results obtained point to an improved accuracy in all metrics tested when dual‐weighing choice of snapshots is combined with POD adaptivity of the trust‐region type. Results of ad‐hoc adaptivity of the POD 4‐D VAR turn out to yield less accurate results than trust‐region POD when compared with high‐fidelity model. Directions of future research are finally outlined. Copyright © 2009 John Wiley & Sons, Ltd.     

A dual‐weighted trust‐region adaptive POD 4‐D Var applied to a finite‐volume shallow water equations model on the sphere

In this paper we study solutions of an inverse problem for a global shallow water model controlling its initial conditions specified from the 40‐yr ECMWF Re‐analysis (ERA‐40) data sets, in the presence of full or incomplete observations being assimilated in a time interval (window of assimilation) with or without background error covariance terms. As an extension of the work by Chen et al. (Int. J. Numer. Meth. Fluids 2009), we attempt to obtain a reduced order model of the above inverse problem, based on proper orthogonal decomposition (POD), referred to as POD 4D‐Var for a finite volume global shallow water equation model based on the Lin–Rood flux‐form semi‐Lagrangian semi‐implicit time integration scheme. Different approaches of POD implementation for the reduced inverse problem are compared, including a dual‐weighted method for snapshot selection coupled with a trust‐region POD adaptivity approach. Numerical results with various observational densities and background error covariance operator are also presented. The POD 4‐D Var model results combined with the trust‐region adaptivity exhibit similarity in terms of various error metrics to the full 4D Var results, but are obtained using a significantly lesser number of minimization iterations and require lesser CPU time. Based on our previous and current work, we conclude that POD 4‐D Var certainly warrants further studies, with promising potential of its extension to operational 3‐D numerical weather prediction models. Copyright © 2011 John Wiley & Sons, Ltd.     

Evaluation of well‐balanced bore‐capturing schemes for 2D wetting and drying processes

We consider numerical solutions of the two‐dimensional non‐linear shallow water equations with a bed slope source term. These equations are well‐suited for the study of many geophysical phenomena, including coastal engineering where wetting and drying processes are commonly observed. To accurately describe the evolution of moving shorelines over strongly varying topography, we first investigate two well‐balanced methods of Godunov‐type, relying on the resolution of non‐homogeneous Riemann problems. But even if these schemes were previously proved to be efficient in many simulations involving occurrences of dry zones, they fail to compute accurately moving shorelines. From this, we investigate a new model, called SURF_WB, especially designed for the simulation of wave transformations over strongly varying topography. This model relies on a recent reconstruction method for the treatment of the bed‐slope source term and is able to handle strong variations of topography and to preserve the steady states at rest. In addition, the use of the recent VFRoe‐ncv Riemann solver leads to a robust treatment of wetting and drying phenomena. An adapted ‘second order’ reconstruction generates accurate bore‐capturing abilities.This scheme is validated against several analytical solutions, involving varying topography, time dependent moving shorelines and convergences toward steady states. This model should have an impact in the prediction of 2D moving shorelines over strongly irregular topography. Copyright © 2006 John Wiley & Sons, Ltd.     

A general approximate‐state Riemann solver for hyperbolic systems of conservation laws with source terms

An approximate‐state Riemann solver for the solution of hyperbolic systems of conservation laws with source terms is proposed. The formulation is developed under the assumption that the solution is made of rarefaction waves. The solution is determined using the Riemann invariants expressed as functions of the components of the flux vector. This allows the flux vector to be computed directly at the interfaces between the computational cells. The contribution of the source term is taken into account in the governing equations for the Riemann invariants. An application to the water hammer equations and the shallow water equations shows that an appropriate expression of the pressure force at the interface allows the balance with the source terms to be preserved, thus ensuring consistency with the equations to be solved as well as a correct computation of steady‐state flow configurations. Owing to the particular structure of the variable and flux vectors, the expressions of the fluxes are shown to coincide partly with those given by the HLL/HLLC solver. Computational examples show that the approximate‐state solver yields more accurate solutions than the HLL solver in the presence of discontinuous solutions and arbitrary geometries. Copyright © 2006 John Wiley & Sons, Ltd.     

A finite volume solver for 1D shallow‐water equations applied to an actual river

This paper describes the numerical solution of the 1D shallow‐water equations by a finite volume scheme based on the Roe solver. In the first part, the 1D shallow‐water equations are presented. These equations model the free‐surface flows in a river. This set of equations is widely used for applications: dam‐break waves, reservoir emptying, flooding, etc. The main feature of these equations is the presence of a non‐conservative term in the momentum equation in the case of an actual river. In order to apply schemes well adapted to conservative equations, this term is split in two terms: a conservative one which is kept on the left‐hand side of the equation of momentum and the non‐conservative part is introduced as a source term on the right‐hand side. In the second section, we describe the scheme based on a Roe Solver for the homogeneous problem. Next, the numerical treatment of the source term which is the essential point of the numerical modelisation is described. The source term is split in two components: one is upwinded and the other is treated according to a centred discretization. By using this method for the discretization of the source term, one gets the right behaviour for steady flow. Finally, in the last part, the problem of validation is tackled. Most of the numerical tests have been defined for a working group about dam‐break wave simulation. A real dam‐break wave simulation will be shown. Copyright © 2002 John Wiley & Sons, Ltd.     

A semi‐implicit numerical method for the free‐surface Navier–Stokes equations

Semi‐implicit methods are known for being the basis of simple, efficient, accurate, and stable numerical algorithms for simulating a large variety of geophysical free‐surface flows. Geophysical flows are typically characterized by having a small vertical scale as compared with their horizontal extents. Hence, the hydrostatic approximation often applies, and the free surface can be conveniently represented by a single‐valued function of the horizontal coordinates. In the present investigation, semi‐implicit methods are extended to complex free‐surface flows that are governed by the full incompressible Navier–Stokes equations and are delimited by solid boundaries and arbitrarily shaped free‐surfaces. The primary dependent variables are the velocity components and the pressure. Finite difference equations for momentum, and a finite volume discretization for continuity, are derived in such a fashion that, after simple manipulation, the resulting pressure equation yields a well‐posed piecewise linear system from which both the pressure and the fluid volume within each computational cell are naturally derived. This system is efficiently solved by a nested Newton type iterative scheme, and the resulting fluid volumes are assured to be nonnegative and bounded from above by the available cell volumes. The time step size is not restricted by stability conditions dictated by surface wave speed, but can be freely chosen just to achieve the desired accuracy. Several examples illustrate the model applicability to a large range of complex free‐surface flows and demonstrate the effectiveness of the proposed algorithm. Copyright © 2013 John Wiley & Sons, Ltd.     

Explicit solution to the exact Riemann problem and application in nonlinear shallow‐water equations

The Riemann solver is the fundamental building block in the Godunov‐type formulation of many nonlinear fluid‐flow problems involving discontinuities. While existing solvers are obtained either iteratively or through approximations of the Riemann problem, this paper reports an explicit analytical solution to the exact Riemann problem. The present approach uses the homotopy analysis method to solve the nonlinear algebraic equations resulting from the Riemann problem. A deformation equation defines a continuous variation from an initial approximation to the exact solution through an embedding parameter. A Taylor series expansion of the exact solution about the embedding parameter provides a series solution in recursive form with the initial approximation as the zeroth‐order term. For the nonlinear shallow‐water equations, a sensitivity analysis shows fast convergence of the series solution and the first three terms provide highly accurate results. The proposed Riemann solver is implemented in an existing finite‐volume model with a Godunov‐type scheme. The model correctly describes the formation of shocks and rarefaction fans for both one and two‐dimensional dam‐break problems, thereby verifying the proposed Riemann solver for general implementation. Copyright © 2007 John Wiley & Sons, Ltd.     

A comparison of the GWCE and mixed P − P1 formulations in finite‐element linearized shallow‐water models

The appearance of spurious pressure modes in early shallow‐water (SW) models has resulted in two common strategies in the finite element (FE) community: using mixed primitive variable and generalized wave continuity equation (GWCE) formulations of the SW equations. One FE scheme in particular, the P ? P₁ pair, combined with the primitive equations may be advantageously compared with the wave equation formulations and both schemes have similar data structures. Our focus here is on comparing these two approaches for a number of measures including stability, accuracy, efficiency, conservation properties, and consistency. The main part of the analysis centres on stability and accuracy results via Fourier‐based dispersion analyses in the context of the linear SW equations. The numerical solutions of test problems are found to be in good agreement with the analytical results. Copyright © 2011 John Wiley & Sons, Ltd.