A Stabilised Nodal Spectral Element Method for
Fully Nonlinear Water Waves
Abstract
We present an arbitraryorder spectral element method for generalpurpose simulation of nonoverturning water waves, described by fully nonlinear potential theory. The method can be viewed as a highorder extension of the classical finite element method proposed by Cai et al (1998) CaiEtAl1998 (), although the numerical implementation differs greatly. Features of the proposed spectral element method include: nodal Lagrange basis functions, a general quadraturefree approach and gradient recovery using global projections. The quartic nonlinear terms present in the Zakharov form of the free surface conditions can cause severe aliasing problems and consequently numerical instability for marginally resolved or very steep waves. We show how the scheme can be stabilised through a combination of overintegration of the Galerkin projections and a mild spectral filtering on a per element basis. This effectively removes any aliasing driven instabilities while retaining the highorder accuracy of the numerical scheme. The additional computational cost of the overintegration is found insignificant compared to the cost of solving the Laplace problem. The model is applied to several benchmark cases in two dimensions. The results confirm the high order accuracy of the model (exponential convergence), and demonstrate the potential for accuracy and speedup. The results of numerical experiments are in excellent agreement with both analytical and experimental results for strongly nonlinear and irregular dispersive wave propagation. The benefit of using a highorder – possibly adapted – spatial discretisation for accurate water wave propagation over long times and distances is particularly attractive for marine hydrodynamics applications.
keywords:
Nonlinear and dispersive free surface waves, Hydrodynamics, Spectral Element Method, Unstructured mesh, Finite Element methods, Highorder discretisation.color \FXRegisterAuthorakapekAK \FXRegisterAuthordbadbDB
1 Introduction
Robust and costefficient timedependent simulation of the propagation and transformation of water waves in both shallow nearshore and deeper offshore areas is a computationally challenging and longstanding scientific problem for ocean, coastal and naval engineering applications. For example, fully nonlinear wave simulations have been subject to research for a long time, and have still not yet entered common coastal and ocean engineering practice. One remaining key challenge is to resolve accurately highly nonlinear and dispersive wave propagation in maritime areas while taking into account varying bathymetry, the geometry of complex structures and their nonlinear interaction with fixed and floating structures. Resolving this problem lead to improved opportunities for using simulations in realistic marine regions as well as enabling experiments in numerical wave tanks of increasing fidelity. Furthermore, it is attractive to design a flexible computational framework that entails use on modern commodity workstations as well as highperformance computing systems. These goals dictate stringent requirements on the design of engineering tools, and this suggests that a generic tool for wave propagation should be based on

a general modelling basis for broadly describing relevant physics,

a generalised computational framework for the numerics, and

software design for efficient mapping to modern and emerging manycore architectures.
Our main objective is to meet each of these requirements via a spectral element based computational framework. Such a framework has already been used for several different applications within marine hydrodynamics such as the fully nonlinear potential flow equations (present work), Boussinesq and shallow water equations EskilssonEngsigKarup2014 (). In this work we focus exclusively on the first two requirements and we pave the way for the fulfilment of the last one, which is being addressed in ongoing work. The use of modern manycore hardware is attractive for acceleration of the highorder spectral element framework and for enabling practical computations of realistic engineering problems Markall2013 (); Markidis01082015 ().
1.1 Choice of modelling basis for description of physics
During the last decades computationally efficient depthintegrated Boussinesqtype models have been widely adopted as essential tools for water wave modelling in the nearshore region; see e.g. MS99 (); Brocchini2013 (). For shorter waves, such as the ones arising in offshore and naval engineering, Boussinesqtype models are not applicable due to the limited accuracy in terms of dispersive and nonlinear properties. For these cases we have to turn to Computational Fluid Dynamics (CFD) models based on the NavierStokes equations, or fully nonlinear potential flow (FNPF) models. The CFD models take viscous effects into account; effects that may be important for breaking waves, load computations, boundary layer effects, etc. Even though CFD is often prohibitively expensive in terms of computational resources when considering simulation of entire sea states EskilssonEtAl2015 (), it is widely used to quantify breaking wave loads and wave runup on structures. CFD models are typically too dissipative as a result of the loworder accuracy imposed by computational limitations for largescale wave simulations. In contrast, already today FNPF models can be used for longtime and largescale wave simulations DucrozetEtAl2007 (); Glimberg2013 (). FNPF solvers can be used for resolution of full sea states in large marine or coastal areas where nonlinear waves interact with fixed or floating structures. The cons of the FNPF models are that they cannot account for nonoverturning breaking waves and viscous effects. For these reasons it can be attractive to combine FNPF models (farfield) with CFD (nearfield) in hybrid modelling approaches for wave structure interaction, cf. PaulsenEtAl2014 (). This hybrid approach enables better simulation of strong nonlinear wave structure interactions in areas where the local wave climates can not be predicted accurately via a FNPF model.
1.2 On the quest towards developing numerical strategies for realworld applications
A review of existing conventional discretisation methods and applications reveals that historically the main emphasis has been on Finite Difference Methods (FDM), boundary element methods (BEM) and finite element methods (FEM) KimEtAll1999 (); HarrisEtAl2014 (). These methods have been designed for the concept of Numerical Wave Tanks (NWT). The main computational bottleneck in all such numerical solvers is the solution of a large linear system. In FDM and FEM the discretisation procedure leads to sparse linear systems due to the local support of discrete operators, while in BEM it is only the domain boundary that needs representation. The discretisation procedure for BEM is based on a surface integral formulation together with Green’s identities. This leads to dense nonsymmetric matrix operators that cannot be solved in a straightforward way with linear asymptotic scaling. There has been some recent progress in bringing the asymptotic cost (scaling rate) down for BEM HarrisEtAl2014 () for both matrixvector multiplications and storage requirements using both highorder basis functions and the Fast Multipole Method (FMM) GreengardRokhlin1987 (). While this strategy can asymptotically achieve linear complexity n ( number of computational nodes) in work effort for the spatial solver, it has a large constant in front of this asymptotic scaling term due to the need of solving a dense linear system of equations. This leaves BEM less efficient compared to volumebased discretisation methods such as FDM and FEM solvers as suggested in WuTaylor1995 (); ShaoEtAl2014 (). We note, that BEM is particularly attractive as a nearfield solver for cases where waves interact with complex geometries ZhouEtAl2015 () and may be combined with a farfield solver such as FEM WuTaylor2003 (). The overall efficiency and scalability of BEM HarrisEtAl2014b () can be compared to efficient and massively parallel free surface hydrodynamics solvers such as EngsigKarupEtAl2011 (); EGNL13 (); ShaoEtAl2014 () which can achieve very high efficiency and scalability using multigridtype methods LiFleming1997 (); EngsigKarup2014 () for arbitrary sized discrete problems, in particular when the (possibly curvilinear multiblock) meshes are logically structured, e.g. as in Glimberg2013 ().
1.3 Stateoftheart in finite element methods for fully nonlinear water waves
Reviews on the stateofthe art of numerical models for freely propagating water waves are given in KimEtAll1999 (); Tanizawa2000 (); DiasBridges2006 (); SpinnekenEtAll2012 (); NimmalaEtAl2012 (). Our scope in the present work is restricted to FNPF solvers and FEM.
The use of FEM for fully nonlinear water waves was pioneered by Wu & Eatock Taylor (1994) WuTaylor1994 (). Since then, the majority of solvers for fully nonlinear potential flow equations have been limited to second (low) order FEM schemes WuTaylor1995 (); GreavesBorthwickTaylor1997 (); MaWuTaylor2001 (); MaWuTaylor2001b () based on a Mixed Eulerian Lagrangian method LonguetHigginsCokelet1976 () for updating the free surface variables. This approach requires expensive mesh update techniques and may suffer from stability problems for deformed meshes RobertsonSherwin1999 (). However, it is particularly well suited for dealing with the interaction between waves and fixed or freely floating bodies WangWu2011 (). To overcome this expensive remeshing problem a less expensive Quasi Arbitrary Lagrangian Eulerian finite element method (QALEFEM) was developed MaYan2006 () and a novel mesh update technique for ALE is described in BouffanaisDeville2006 ().
Very little work has been done based on purely Eulerian formulations, e.g. the classical approach where a transformed formulation is used. Previous studies based on using the transformation CaiEtAl1998 (); WesthuisAndonowati1998 (); ClaussSteinhagen:1999tr (); TurnbullEtAl2003 () have been limited to two spatial dimensions and flat bathymetries for wavestructure interaction applications. These are essentially limited to the inclusion of bottommounted possibly surfacepiercing structures, despite that the method may be used as an efficient base solver for very large domains and from shallow to deep waters. More details about free surface solvers based on FEM can be found in Zienkiewicz2014 (); TurnbullEtAl2003 (); WangWu2011 ().
While there are very few studies on highorder finite element methods for water wave applications, it is well known that large efficiency gains for timedependent wave problems can be achieved by the use of highorder accurate methods KO72 (). One attractive class of methods are the Spectral Element Methods (SEM) EESHB05 (); DFM02 (); KS2005 () which rely on the strong theoretical foundations of Spectral Methods CHQZ2006 (); HGG07 (). The Spectral Element Method was first used by Patera (1984) PAT84 () for fluid dynamics problems and has since then gained popularity. It combines the best properties of spectral methods and classical finite element methods, namely, high accuracy and flexibility in the spatial representation of domains. For smooth problems, the use of highorder discretisation is generally an efficient way to balance accuracy and cost, since a highorder discretisation allow for more coarse spatial representation compared to a loworder method. Highorder methods may also reduce the robustness of explicit timestepping methods due to global restriction on the stable time step sizes, although, in recent works EGNL13 (); EskilssonEngsigKarup2014 () it is demonstrated that for certain wave models that require operator inversion in the timestepping, highorder discretization combined with explicit time stepping methods need not reduce robustness  not even for unstructured mesh methods – and thereby provide a strong basis for efficient tools. Furthermore, with the geometric flexibility provided by adaptive meshes and the support for costefficient simulation through  and refinement strategies, the SEM is suitable for largescale computing due to accurate temporal integration over long times via highorder discretisation that lead to small dispersive and dissipative errors. However, a key challenge for unstructured solvers such as FEM and SEM is achieving work effort that is linearly proportional to the number of computational processing nodes (weak scalability) for unstructured meshes used to represent general domains. This requires the design of efficient and scalable preconditioned iterative methods CaiEtAl1998 () and of data structures that map efficiently to modern manycore architectures GoddekeEtAl2008 ().
1.4 Paper contributions
The main objective of this work is the proper design and validation of the SEM framework for – ultimately largescale – dispersive and nonlinear water wave propagation. This effort is crucial to enable engineering analysis of water waves in marine regions. This paper proposes for the first time a stabilised highorder spectral element method for solving the fully nonlinear potential flow equations. A rigorous assessment of the stabilised numerical model is carried out using numerical experiments and known benchmarks in two space dimensions.
2 Governing equations
For the description of inviscid and irrotational fluid flows we introduce the set of fully nonlinear and dispersive free surface equations described by a twovariable potential flow formulation. This formulation can be derived from the NavierStokes equations, cf. EGNL13 ().
Let the fluid domain () be a bounded, connected domain with piecewise smooth boundary and introduce restrictions to the free surface and the bathymetry . Let be the time domain. The mathematical problem is to find a scalar velocity potential function and to determine the evolution of the free surface elevation . The notations are illustrated in Figure 1.
The Eulerian description of the unsteady kinematic and dynamic free surface boundary conditions can be expressed in the Zakharov form ZAK1968 (). Find such that
(1a)  
(1b) 
We have introduced the ’’ symbol to denote functionals defined only on the free surface plane. The vertical component of the velocity is determined by first solving a Laplace problem
(2a)  
(2b)  
(2c) 
where describes the still water depth.
We are interested in a set of governing equations that can be used as a basis for efficient simulations with support for representing structures accurately. A basis for efficient simulations is the classical transformation of the vertical coordinate
(3) 
where is the height of the water column above the bottom. This transforms the fluid domain to a timeindependent computational domain at the expense of timevarying coefficients. The notations are illustrated in Figure 2.
The main drawback of using the transformation is that it needs to be nonsingular, and thus exclude the geometric modelling of breaking waves. For general wavestructure problems, this restriction can be removed by discretising and solving the Laplace problem (2) directly. Following CaiEtAl1998 (), we express the transformed system in a form where variable depth is accounted for. Let be the timeindependent computational domain . The Jacobian of the map is then
(4a)  
enabling the transformed system to be expressed in the differential form  
(4b)  
where is introduced and the symmetric coefficient matrix is  
(4c) 
The artificial scalar velocity function contains all information about the flow kinematics in the entire fluid volume. The velocity field can be determined from using the relation . All metric coefficients in (4a) can be evaluated from the known twodimensional free surface and bottom positions at given instants of time. It is possible to discretise (4b) after completing the differentiations, cf. EBL08 (); WesthuisAndonowati1998 (), however, this increases the complexity of this formulation.
2.1 Boundary conditions
For the solution of the Laplace problem at every time step, the following free surface boundary condition is specified
(5) 
while at vertical boundaries, impermeable wall boundary conditions are assumed
(6) 
where is an outward pointing unit normal vector at . Due to symmetry at wall boundaries the domain boundary conditions at free surface variables imposed are
(7) 
2.2 Wave generation and absorption zones
We employ a generalpurpose embedded penalty forcing technique EGNL13 () similar to the technique described in Cointe1989 (). This technique can be used for both wave generation of regular and irregular waves as well as absorption in a numerical wave tank setup. It can be derived from the relaxation method described in LD83 () to avoid the preprocessing step and turn it into an equivalent source term to the governing equations
(8) 
where is the nonlinear function for the PDE, a free surface state variable (, in Eq. (1)), the relaxation functions can be defined to avoid minimal reflections ENG06 (), and the source function defines the analytical representation of the wave input signal to be generated. We choose equal to time step size .
3 Numerical discretisation
The governing partial differential equations are discretised in a generic computational framework based on the method of lines, where first a semidiscrete system of ordinary differential equations is formed by spatial discretisation using a nodal Spectral Element Method. We present the 2D formulation before presenting results of 2D numerical experiments.
3.1 Weak Galerkin formulation and discretisation
We form a partition of the domain to obtain a tessellation of consisting of nonoverlapping shaperegular elements such that with denoting the ’th element. We introduce the finite element approximation space of continuous, piecewise polynomial functions of degree at most , .
3.1.1 Unsteady free surface equations
The weak formulation of the free surface equations takes the following form. Find where such that
(9a)  
(9b) 
for all . We introduce the finitedimensional approximations
(10) 
where is the set of global finite element basis functions with cardinal property at mesh nodes with the Kronecker Symbol. Substitute these expressions into (9) and choose . The discretisation in one spatial dimension becomes
(11a)  
(11b) 
where the following global matrices have been introduced
(12) 
The gradients of the free surface state variables are recovered as described in Section 3.1.3. Temporal integration of (11) is performed using an explicit fourthorder RungeKutta method.
3.1.2 Spatial discretisation of the Laplace problem
Consider the discretisation of the governing equations for the transformed Laplace problem (4b). We seek to construct a linear system of the form
(13) 
where is the total degrees of freedom in the discretisation.
The starting point is the weak formulation of the symmetric formulation (4b) that can be expressed as: find such that
(14) 
for all where the boundary integrals vanish at domain boundaries where impermeable walls are assumed. The discrete system operator is defined via domain decomposition as
(15) 
The elemental integrals are approximated through the change of variable
(16) 
where is the Jacobian of the affine mapping and is the computational reference element. The global assembly of this operator preserves symmetry, and the resulting linear system is then modified to impose the Dirichlet boundary conditions (5) at the free surface. Typically, is a large and sparse operator with a narrow band structure of nonzero elements in two space dimensions obtained by a proper permutation of the rows and numbering of nodes to reduce fillin in the factorisation of Gaussian Elimination procedures. By using a symmetric reverse CuthillMackee permutation, the bandwidth of the sparse matrix is minimised and the system can be efficiently solved by a sparse direct Gaussian elimination procedure. In 2D this leads to optimal and scalable work effort and is used in this work. For sparse symmetric positive definite systems such as Eq. (13), the iterative preconditioned conjugate gradient solver is an attractive choice when system sizes become large CaiEtAl1998 () as would be expected for typical applications in three space dimensions, since convergence is guaranteed and memory footprint is minimal.
3.1.3 A generic technique for gradient recovery
The gradients of the globally piecewise continuous basis functions will be discontinuous across element interfaces in the classical sense. To guarantee global continuity of derivatives a gradient recovery procedure can be used. Several gradient recovery techniques are reviewed in HintonCampbell1974 (); HawkenTownsendWebster1991 (); SriramEtAl2010 (). In this work, a global gradient recovery technique is used within our spectral element framework.
The global approximation of components of the horizontal first derivative as a function is expressed as
(17) 
By a global Galerkin projection
(18) 
we can generate the linear systems of equations
(19) 
to recover the coefficients of the expansion (17). This procedure is similar to the one described by HawkenTownsendWebster1991 () and used in other models, cf. WuTaylor1994 (); RobertsonSherwin1999 (); TurnbullEtAl2003 ().
Remark: In the FEM model described in Westhuis2001 () the gradient recovery step is related to stability of the numerics. It was found that a global projection method may lead to instability when using FEM. As we shall see in the numerical experiments in Section 5 we do not reach the same empirical conclusion.
Using this gradient recovering technique and after the solution of (2), it is possible to estimate the vertical free surface velocity to obtain closure in the free surface problem (1b). The weak formulation is
(20) 
From this we can construct a linear system of the form
(21) 
where the discrete operators takes the form
(22) 
From the solution of (21) we obtain the vertical free surface velocities
(23) 
where denotes an index set consisting of global numbers for the free surface nodes.
3.1.4 Nodal prismatic Lagrange finite elements in two space dimensions
We consider first construction of elements in two space dimensions, needed for representing the solutions of the governing equations. The expansions can be based on quadrilateral elements in the timeconstant computational domain. These elements are formed by a subdivision of the horizontal free surface plane with possible irregular sized nonoverlapping elements. Each of these elements can then be extended in the vertical from the surface to the bottom to form the quadrilateral elements. A spectrally accurate multivariate hierarchical polynomial expansion can be constructed by a tensor product of onedimensional basis functions. This results in a quadrilateral element . We introduce the element basis functions
(24) 
where is the ’th order orthonormal Jacobi polynomial on the interval with orthogonality property
(25) 
These polynomials can be evaluated efficiently using a simple recurrence relation HGG07 (); Kopriva2009 (). The basis (24) have arbitrary polynomial orders in the horizontal plane and in the vertical plane. This makes it possible to tune the orders of the approximations to balance accuracy and efficiency needs in simulations.
For interpolating polynomials, the Unisolvence Theorem guarantees a unique connection between the hierarchical polynomial (modal) expansion and the corresponding Lagrange polynomial (nodal) expansion. Thus, for all , , and , and for each element , we have
(26) 
which defines a relationship between modal and nodal coefficients in the form HW08 ()
(27) 
where a nonsingular generalised Vandermonde matrix has been introduced. This connection can be exploited together with the duality in polynomial representation to define the local Lagrange basis functions with property and their derivatives
(28) 
and the matrices defined in terms of the first derivatives of the modal basis functions
(29) 
For use with integration on the local elements, the local nodal mass matrix is introduced
(30) 
where orthonormality of the basis functions is exploited to avoid the use of discrete quadrature rules in the constructions, leaving the implementations quadraturefree.
By evaluating expressions such as (28) at the chosen inter element node distribution we obtain the generic arbitraryorder elemental operators
(31) 
that can be used in connection with the chain rule to define derivatives in physical space
(32) 
with discrete counterparts based on collocation expressed in the form
(33) 
The choice of the nodal distribution on the simplexes in two horizontal dimensions can be based on explicit construction of the nodal distribution set using, e.g., a blend and warp procedure War06 (). This is used to define the highorder nodal basis functions WarburtonSherwinKarniadakis1999 ().
3.2 Generalised local element matrices via quadraturefree matrixbased operations
All global operators can be assembled from generalised local elemental operators. Consider global integrals in the general form
(34) 
Approximate each of the integrands using finite element basis functions of the form
(35) 
The global integrals can be reduced to local integrals through domain decomposition
(36) 
and by inserting the approximated integrands, we obtain local integrals of the form
(37) 
Expand the approximate integrands using nodal expansions such that
(38) 
In the special cases, where and are differential operators, e.g.
(39) 
the local integrals in (38) can be evaluated by exploiting that for nodal differentiation matrices , cf. CHQZ2006 (). An implementation of the local elemental operator following the derivation just outlined, enables generic implementation in a single framework that has support for adaptivity to balance accuracy and cost in computations. Furthermore, the integration of the integrands approximated using polynomial basis functions are without quadrature errors when quadraturefree matrixbased operations are carried out using evaluations on sufficiently fine meshes as described next.
3.3 Removal of quadrature errors via higherorder numerical integration
Numerical integration of nonlinear terms is handled via supercollocation KK03 () and is expressed in the general form
(40) 
where we have introduced the interpolation operator that maps the representation of a function on coarse mesh (subscript ) to that of a fine mesh (subscript ). Exploiting the result in (38), the integration (projection) of the polynomial representation can be done without quadrature errors via mass matrix based operations for integration of the form
(41) 
where and . The computational cost of this integration step has complexities and due to the involved matrix operation. The relation between size of the coarse and fine variables is for exact integration of the polynomial representation of the coarse space basis functions when quartic terms are present. Therefore, the increase in cost is about four times for these operations.
3.4 Dealing with aliasing by spectral filtering
The solution of governing equations that contain strongly nonlinear terms may pose a challenge for maintaining both accuracy and numerical stability in time when simulations are marginally resolved. In such cases, numerical instability may be related to aliasing effects. Aliasing results from evaluating interpolated products of functions, which when represented with insufficient or marginal resolution introduces errors in the functional representation. Possible remedies KK03 (); Mengaldo201556 () focus on increasing resolution or introducing a proper stabilisation strategy to artificially dissipate such errors. While the former is simple, it is generally not considered feasible for largescale systems or for longtime integration.
In this work, we employ a spectral filtering strategy exploiting the duality in the local element representation, cf. (26). On the ’th element the filtered local solution can be expressed as (assuming an order expansion in one space dimension here)
(42) 
The filtering is applied only for the timedependent free surface variables and . An exponential filter HK08 (); HW08 () can be used with cutoff lowpass filter index , to only affect the highest modes , such that
(43) 
For example, by choosing the parameters a mild damping is achieved which gently removes five percent of the energy from only the highest mode in the basis. The filtering is done on a per element basis using a matrixvector product
(44) 
which has a work complexity of with total work effort proportional to . This local filtering matrix is constructed and used on all elements and applied repeatedly as necessary. We use the model basis in one space dimension KS99 ()
(45) 
to avoid introducing interface jumps or affect the mean through only filtering the higher modes (). If the filter is applied gently, i.e. by removing very little ’energy’ from the highest modes only, spectral accuracy can be recovered. Excessive filtering may reduce the convergence rate of the method albeit improve stability.
4 Numerical properties of the model
4.1 Temporal stability
For general schemes, the connection between nonlinear instability and explicit timeintegration is important to understand DiasBridges2006 (). To solve the governing equations efficiently, an explicit time integration method is preferred. Explicit schemes come with conditional stability in the form of the global CFL condition that dictates an upper bound for the choice of stable time step sizes where the constant . It is wellknown that a main challenge for many highorder Spectral Element Methods is an unattractive scaling of the modulus of the eigenvalues of the discrete operators of the form where is expansion order and is the highest order of differentiation operator in the evolution equations. Typically, the constant is dependent on the minimum mesh size for an element in the mesh. This property may pose a severe problem in the accurate local representation of geometric features with small elements.
As discussed and shown in EGNL13 () for a linear finite difference scheme of the same governing equations, the stable time step size for explicit schemes has an upper bound given by the CFL condition , where is here still water depth and the gravitational acceleration constant. This bound is only dependent on the scale of the physics (still water depth) and the resolution chosen in the vertical. This result is demonstrated and shown in Figure 3, following EBL08 () via discretisation and numerical eigenvalue analysis of the linearised system. This property is similar to the property inherent in several other wave models as described in the works RobertsonSherwin1999 (); ENG06 (); EHBM06 (); EHBW08 (); EGNL13 (); EskilssonEngsigKarup2014 (). This property implies that the CFL condition is tamed WarburtonHagstrom2008 () in the sense that the time step size is not dictated by the numerics but only the physics (depth). Along the same line of the experiments described in EGNL13 (), we find that there are only small changes in the CFL properties for nonlinear problems and hence small elements in the mesh do not impose any severe time step size restriction.
4.2 Linear accuracy and dispersion properties
The accuracy of the numerical model is compared to the theoretical solution to the system of equations that arise when the system is subject to the assumption of smallamplitude waves (assume where wave height relative to wave length ). The theoretical solution for linear progressive monochromatic waves in one space dimension is given by SJ01 ()
(46) 
with linear dispersion relation from Stokes theory , where is the wave number, is the gravitational acceleration (assumed as 9.81 ) and is the angular velocity with being the wave period.
In Figure 4 (a) results of  and refinement strategies are presented and highlight an important advantage of the SEM. When solutions are smooth, it is possible to use highorder basis functions to improve the costefficiency of the method by allowing fewer degrees of freedom to be used to achieve the same level of accuracy of lower order methods. Most previous works have focused exclusively on classical FEM methods based on piecewise linear approximations, which has a convergence rate that matches the curve for linear () basis functions.
The decisive criterion for choosing between different numerical strategies is to understand what is the amount of work (cost) to achieve a given level of accuracy. This is illustrated via results obtained via a semioptimal bruteforce procedure in Figures 4 (b) and (c) where the work effort has been minimised^{1}^{1}1Tests done using a sequential proofofconcept code on a laptop with a 2,3 GHz Intel Core I7 processor and 16 GB 1600 MHz DDR3 RAM. with respect to time to be as numerically efficient for a fixed accuracy level of 1 error in surface elevation. With no significant dissipation in the scheme, this implies that the error measure mainly numerical dispersion error. The results in Figure 4 (d) show that for our current sequential proofofconcept implementation there is speedup of approximately x2 for short time and close to x6 for longer times, achievable by switching from a ()scheme to a scheme. These gains can be improved further to close to x3 for short time and close to x37 for longer times, achievable by switching from a ()scheme to a ()scheme. As expected (cf. KO72 ()), the longer a simulation the more the possible gain in terms of numerical efficiency, and for even higher orders, there are additional gains although they end up being marginal for the modest accuracy requirement chosen.
Clustering mesh nodes more densely closer to the free surface can improve accuracy in linear dispersion without increasing the CPU time significantly MaWuTaylor2001 (); BinghamZhang2007 (); EBL08 (). This is confirmed by results presented in Figure 5. The results show that the use of a highorder method or clustered vertical distribution of loworder elements is a must for an accurate approximation of dispersion in applications where is large (short wave length relative to depth). This highlights that one can tune the accuracy by proper choice of discretisation parameters. An important implication of these results is that the vertical node distribution can be used to control the range of validity of the model in terms of dispersive properties, i.e. a numerical truncation counterpart to the analytic truncation used in Boussinesqtype models. Interestingly, flexibleorder finite difference solvers with cosineclustered vertical distributions of nodes appears to be superior in linear dispersion compared to single layer SEM for equal number of vertical nodes. However, multilayer SEM with cosine clustered elements of vertical order almost match the FDM, cf. (EGNL13, , Figure 2.6, page 59).
4.3 On nonlinear accuracy, stability and kinematics properties
The accurate computation of kinematics is essential for load predictions in wavestructure applications, e.g. for offshore foundations of wind turbines. For nonlinear waves, exact stream function (SF) wave solutions of permanent form Dean1965 () based on assuming a flat sea bed can be used to assess the accuracy with respect to variation in dimensionless nonlinearity () and dispersion () parameters. This is done by solving the Laplace problem first. Then, from the scalar velocity potential solution we can calculate the vertical free surface velocity and compare with exact results. Numerical results for linear, weakly nonlinear and strongly nonlinear SF waves in combinations of shallow, intermediate and deep waters are presented in Figures 6 and 7. Here we use only one layer of elements in the vertical and a fixed number of elements in the horizontal direction. Convergence is achieved by the variation of the polynomial order to achieve fast convergence. All tests show convergence with increasing resolution as expected. When depth or nonlinearity increases more resolution is required. Similar tests were carried out for a flexibleorder finite difference model in BinghamZhang2007 (). An immediate conclusion is that the SEM has similar resolution requirements as the corresponding finite difference solver to match the order of accuracy for nonlinear applications. This highlights that from a purely algorithmic (implementation independent) viewpoint there seems to be no significant tradeoffs in introducing geometric flexibility through this choice of discretisation. In Figure 8, the accuracy of the kinematics computations is shown for intermediate depth and deep water for very steep nonlinear stream function waves. Excellent agreement is found between exact and computed results for both intermediate and deep waters, which is difficult to represent in conventional wave propagation models due to lack of resolution.
5 Numerical experiments
We examine different test cases and benchmarks, that inspect different properties of the numerical model that serves as validation.
5.1 Stabilised nonlinear wave propagation of stream function waves
To test the robustness and accuracy of the numerical method, we consider the strenuous test case of propagating nonlinear stream function waves near the theoretical limit of wave steepness and nonlinearity. Few numerical schemes can resolve such waves accurately. Results from representative numerical experiments are presented in Table 1. The table shows how a standard Galerkin scheme unsurprisingly fails to be stable for all temporal resolutions chosen. Improvement in stability is achieved by using higher order numerical integration (cf. Section 3.3) integrands in the Galerkin integrals. This is done by interpolating the integrands to a finer mesh on the reference element corresponding to a basis with polynomial basis expansion order of resulting in terms of less than quartic order nonlinear to be subject to overintegration. This implies that exact integration is used for all nonlinear terms to within the order of accuracy of the numerical scheme. This is compared to using no overintegration but instead trying to stabilise only using a gentle spectral filtering strategy that caps off the highest 5% in the highest modes of the modal expansion. Both of these strategies fail to be stable. Instead, if the two strategies are combined, where the overintegration effectively removes any aliasing in the evaluation of the strongly nonlinear terms, the mild spectral filter dissipates just enough energy for the model to stabilise completely. We find that overintegration is only needed in the freesurface equations, leading to a marginal increase in computational cost, which is anyway driven by the Laplace solver. This is highlighted in the table, where the time stepping cost is only increased by approximately 15% when using overintegration and spectral filtering in comparison with a standard Galerkin formulation based on expansion order . This additional cost of overintegration would be even less significant for larger simulations.
Accuracy  Cost/  

40  80  160  
No filter +  NaN  NaN  NaN  1.00  
No overintegration  NaN  NaN  NaN  
No filter +  1.3608e03  7.5021e04  7.6731e04  1.02  
Overintegration  NaN  NaN  NaN  
Filter +  NaN  5.6019e04  8.0508e03  1.12  
No overintegration  NaN  NaN  1.2705e03  
NaN  NaN  7.8374e03  
NaN  NaN  NaN  
Filter +  1.3943e03  7.0651e04  1.0102e03  1.15  
Overintegration  7.4032e03  4.3313e03  7.0332e03  
7.2826e02  5.7642e02  7.5093e02 
5.2 Convergence tests and highorder accuracy
To confirm the highorder accuracy of the model and evaluate the influence of spectral filtering we have carried out convergence tests using the exact nonlinear stream function wave solutions for parameters (dispersion) and , and of maximum wave steepness (nonlinearity). The results for tests of the proposed Galerkin scheme with overintegration without spectral filtering and with spectral filtering using a cap of 1% of the highest modal mode are presented in Figure 9. The results confirm the highorder convergence for the spatial spectral element discretisation. Particularly, this is clear for the mildly nonlinear wave. With increasing nonlinearity more spatial resolution is required to accurately resolve all harmonic modes of the solution. The gentle filtering is found to reduce accuracy and has some detrimental effect on the convergence. For increasing order of the local basis functions, these effects become less significant.
5.3 Harmonic generation over a submerged bar
We consider the classical benchmark for wave transformation due to a submerged bar to test the accuracy of the SEM model. The test is often used for validation of deterministic dispersive and nonlinear wave models since it can be compared to experimental laboratory data and other known results in the literature, e.g., see GK99 (); MS99 (). The numerical wave tank is illustrated in Fig. 10.
The experiment was originally proposed in BB94 () and subsequently an equivalent scaled experiment was carried out and described in LKK94 (). We consider the setup for Case A in the original experiment. The weakly nonlinear input wave is generated in the numerical model using a regular stream function solution at undisturbed depth m with a wave height of cm and a wave period of s. The input wave is generated and propagates towards the submerged bar on a constant bottom. During the propagation over the bar, the wave will undergo a transformation resulting in a steepening and shortening of wavelengths due to nonlinear shoaling effects. At the top of the bar, the bound harmonics will be released as free harmonics (harmonic generation) decomposing the wave into shorter waves that propagate freely. Thus, to attain high accuracy in the calculations we need to use a model that can handle nonlinear wavewave interactions and have accurate dispersion properties to capture the correct wave speed of the free harmonics after the bar. Taking advantage of the unstructured SEM the elements are of similar sizes but adjusted to have the interfaces positioned where the depth function has kinks in the first order gradients. Thus the transformation is also local to the elements and entirely valid throughout the domain. The results presented in Figure 11 are based on a deterministic simulation where we have used 103 elements in the horizontal and one vertical layer with a multivariate basis of polynomial order in both the horizontal and vertical directions. A CFL condition with Courant number is used for defining the time step size. All results are found to be in excellent agreement with the experimental data with some qualitatively minor differences in phase between experimental and computed results which compare well to other published results.
5.4 Irregular waves shoaling on a slope
Mase and Kirby MaseKirby1992 () conducted laboratory experiments of shoaling and breaking of irregular waves on a constant slope. The water depth offshore the slope is 0.47 m and the slope is 1/20. The incoming irregular waves were generated by a PiersonMoskowitz (PM) spectrum with peak frequency of 1.0 Hz. This set of experiments has been used for testing Boussinesq models in e.g. WK95 (); Shietal12 ().
As the SEM model presently is limited to nonbreaking waves we follow the approach of WK95 (), i.e. the constant slope part of the wave tank is truncated and exchanged with a flat bottom at a depth such that no breaking occurs. For the SEM model we limit the still water depth to be no less than 0.19 m. The numerical wave tank is illustrated in Fig. 12. The wave gauges are located at and m, respectively. Wave generation and absorption is performed with 4 m long relaxation zones. Further, as the random phase angle of the spectrum is unknown, the incoming wave train is generated in the numerical model based on a FFT of the measured waves at wave gauge 1 located at the toe of the slope. The wave kinematics is given by superposition of the linear wave modes.
The computational domain is partitioned into elements of length 0.1 m with an expansion basis of polynomial order in the horizontal and vertical directions. Figure 13 shows the simulated and experimental free surface elevation at the different wave gauge locations. There is a reasonably good fit to the experimental data, with some minor differences in peak and through amplitudes. There are two main reasons for the discrepancy. First, the shorter waves of the PM spectra would require a very high polynomial order in order to be properly resolved and are thus diffused. Secondly, the assumption of linear superposition in the generation zone is not correct, giving that the incoming wave train does not exactly match the recorded.
6 Conclusions
We have presented a spectral element model for simulation of fully nonlinear water wave propagation. The main advantages of using the spectral element method is the opportunity for balancing high accuracy with unstructured meshes which can be adapted to geometry of arbitrary shape (sharp corners, curvilinear features, etc.) or features of the solution such as large gradients. The spectral element method’s dual roads towards convergence, namely  and adaptivity, allows for balancing accuracy and cost effectively.
A twodimensional spectral element model for fully nonlinear potential flow is implemented in a general computational framework using quadraturefree construction of local element operators. We have used arbitraryorder multivariate Lagrange (nodal) basis functions in space and an explicit fourth order RungeKutta method in time. The explicit time stepping is efficient as the model has a bounded eigenspectrum, and the stable time step sizes are governed by still water depth and the vertical resolution. The model is stabilised by using overintegration to effectively reduce aliasing errors and mild spectral model filtering to add some artificial viscosity to secure robustness for marginally resolved flows.
The proposed model was shown to have a convergence rate of order , although it is difficult to keep the sharp convergence rate for the most nonlinear waves due to high resolution requirements. Numerical experiments demonstrate that the stabilised model is both robust and accurate. It was shown how the spectral accuracy can be used to substantially reduce the number of degrees of freedom per wavelength and compare well to properties that can also be achieved with finite difference time domain schemes, but with the additional benefit of geometric flexibility. Also, we illustrate how the order of the basis function in vertical dimension can be used to control the range of validity of the model in terms of dispersive properties, i.e. a numerical truncation counterpart to the analytic truncation used in standard Boussinesqtype models. While the methodology is efficient in two space dimensions, particular attention is to be given to further improve numerical efficiency via efficient preconditioning methods that maintain high efficiency for general unstructured grids.
In ongoing work, we aim at considering advanced nonlinear hydrodynamics problems by extending the current framework to also handle moving and floating objects. The present model needs to be further improved to handle runup for calculations in the swash zone and a breaking wave model needs to be included for realistic applications. The proposed methodology in two space dimensions can be extended without conceptual modifications for three space dimensions, paving the road for marine hydrodynamics applications in 3D for offshore structures. Main challenges for enabling efficient computations in three space dimensions would be to (i) switch from a direct solver to an iterative solver strategy and identify an efficient and scalable preconditioning strategy, and (ii) via proper software design enable efficient mapping to modern and emerging manycore architectures for accelerated performance of the computational framework.
References
References
 [1] S. Beji and J. A. Battjes. Numerical simulation of nonlinearwave propagation over a bar. Coastal Engineering, 23:1–16, 1994.
 [2] H. B. Bingham and H. Zhang. On the accuracy of finitedifference solutions for nonlinear water waves. J. Engineering Math., 58:211–228, 2007.
 [3] E. Bouffanais and M. O. Deville. Mesh update techniques for freesurface flow solvers using spectral element method. J. Sci. Comput., 27(13):137–149, June 2006.
 [4] M. Brocchini. A reasoned overview on boussinesqtype models: the interplay between physics, mathematics and numerics. Proc. Roy. Soc. London Ser. A, 469(2160):1–27, 2013.
 [5] X. Cai, H. P. Langtangen, B. F. Nielsen, and A. Tveito. A finite element method for fully nonlinear water waves. J. Comput. Phys., 143:544–568, July 1998.
 [6] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods  Fundamentals in single domains. Springer, 2006.
 [7] G. Clauss and U. Steinhagen. Numerical simulation of nonlinear transient waves and its validation by laboratory data. International Journal of the International Offshore and Polar Engineering Conference (ISOPE), Brest, France, III:368–375, 1999.
 [8] R. Cointe. Nonlinear simulation of transient free surface flows. In In proceedings of the 5th International Conference in Numerical Ship Hydrodynamics, 1989.
 [9] R. G. Dean. Stream function representation of nonlinear ocean waves. J. Geophys. Res., 70:4561–4572, 1965.
 [10] M. O. Deville, P. F. Fischer, and E. H. Mund. High Order Methods for Incompressible Fluid Flow. Cambridge University Press, 2002.
 [11] F. Dias and T. J. Bridges. The numerical computation of freely propagating timedependent irrotational water waves. Fluid Dynam. Res., 38(12):803–830, 2006.
 [12] G. Ducrozet and P. Ferrant. Rogue waves in largescale fullynonlinear highorderspectral simulations. In Proc. 22nd International Workshop on Water Waves and Floating Bodies (IWWWFB), Croatia, 2007.
 [13] A. P. EngsigKarup. Unstructured Nodal DGFEM solution of highorder Boussinesqtype equations. PhD thesis, Department of Mechanical Engineering, Technical University of Denmark, 2006.
 [14] A. P. EngsigKarup. Analysis of efficient preconditioned defect correction methods for nonlinear water waves. Int. J. Num. Meth. Fluids, 74(10):749–773, 2014.
 [15] A. P. EngsigKarup, H. B. Bingham, and O. Lindberg. An efficient flexibleorder model for 3D nonlinear water waves. J. Comput. Phys., 228:2100–2118, 2009.
 [16] A. P. EngsigKarup, L. S. Glimberg, A. S. Nielsen, and O. Lindberg. Fast hydrodynamics on heterogenous manycore hardware. In Raphäel Couturier, editor, Designing Scientific Applications on GPUs, Lecture notes in computational science and engineering. CRC Press / Taylor & Francis Group, 2013.
 [17] A. P. EngsigKarup, J. S. Hesthaven, H. B. Bingham, and P. Madsen. Nodal DGFEM solutions of highorder Boussinesqtype equations. J. Engineering Math., 56:351–370, 2006.
 [18] A. P. EngsigKarup, J. S. Hesthaven, H. B. Bingham, and T. Warburton. DGFEM solution for nonlinear wavestructure interaction using boussinesqtype equations. Coastal Engineering, 55:197–208, 2008.
 [19] A. P. EngsigKarup, M. G. Madsen, and S. L. Glimberg. A massively parallel GPUaccelerated model for analysis of fully nonlinear free surface waves. Int. J. Num. Meth. Fluids, 70(1), 2011.
 [20] C. Eskilsson and A. P. EngsigKarup. On devising boussinesqtype models with bounded eigenspectra: One horizontal dimension. J. Comput. Phys., 271(0):261 – 280, 2014. Frontiers in Computational Physics Modeling the Earth System.
 [21] C. Eskilsson, A. P. EngsigKarup, S. J. Sherwin, J. S. Hesthaven, and L. Bergdahl. The next step in coastal numerical models: spectral/hp element methods? In Proceedings of the WAVES2005 Conference, Madrid, 2005.
 [22] C. Eskilsson, J. Palm, J. P. Kofoed, and E. FriisMadsen. CFD study of the overtopping discharge of the wave dragon wave energy converter. In Proc. 1st International Conference of Renewable Energies Offshore. ASCE, 2015.
 [23] S. L. Glimberg. Designing Scientific Software for Heterogenous Computing  With application to largescale water wave simulations. PhD thesis, Department of Applied Mathematics and Computer Science, Technical University of Denmark, Kongens Lyngby, Denmark., 2013.
 [24] M. F. Gobbi and J. T. Kirby. Wave evolution over submerged sills: tests of a highorder boussinesq model. Coastal Engineering, 37:57–96, 1999.
 [25] D. Göddeke, R. Strzodka, J. MohdYusof, P. McCormick, H. Wobker, C. Becker, and S. Turek. Using GPUs to improve multigrid solver performance on a cluster. Int. J. Comput. Sci. Eng., 4(1):36–55, 2008.
 [26] D. M. Greaves, G. X. Wu, A. G. L. Borthwick, and R. Eatock Taylor. A moving boundary finite element method for fully nonlinear wave simulations. J. Ship Res., 41(3):181–194, 1997.
 [27] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. Journal of Computational Physics, 135(2):280 – 292, 1997.
 [28] J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. A comparison of methods in fully nonlinear boundary element numerical wave tank development. In 14émes Journées de l’Hydrodynamique, 2014.
 [29] J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. Fast integral equation methods for fully nonlinear water wave modeling. In Proceedings of the Twentyfourth (2014) International Ocean and Polar Engineering Conference, 2014.
 [30] D. M. Hawken, P. Townsend, and M. F. Webster. A comparison of gradient recovery methods in finiteelement calculations. Communications in Applied Numerical Methods, 7(3):195–204, 1991.
 [31] J. S. Hesthaven, S. Gottlieb, and D. Gottlieb. Spectral Methods for TimeDependent Problems. Cambridge Monographs on Applied And Computational Mathematics 21. Cambridge University Press, Cambridge, UK, 2007.
 [32] J. S. Hesthaven and R. M. Kirby. Filtering in legendre spectral methods. Math. Comp., 77(263):1425–1452, 2003.
 [33] J. S. Hesthaven and T. Warburton. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer, 2008.
 [34] E. Hinton and J. S. Campbell. Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Num. Meth. Engng., 8(3):461–480, 1974.
 [35] G. E. Karniadakis and S. J. Sherwin. Spectral/hp element methods for CFD. Oxford University Press, 1999.
 [36] G. E. Karniadakis and S. J. Sherwin. Spectral/hp element methods for computational fluid dynamics. Oxford University Press, 2 edition, 2005.
 [37] C. H. Kim, A. H. Clément, and K. Tanizawa. Recent research and development of numerical wave tanks  a review. Int. J. Offshore and Polar Engng., 9(4):241–256, 1999.
 [38] R. M. Kirby and G. E. Karniadakis. Dealising on nonuniform grids: algorithms and applications. J. Comput. Phys., 191:249–264, 2003.
 [39] D. Kopriva. Implementing Spectral Methods for Partial Differential Equations  Algorithms for Scientists and Engineers. Springer, 2009.
 [40] H.O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24:199–215, 1972.
 [41] J. Larsen and H. Dancy. Open boundaries in short wave simulations  a new approach. Coastal Engineering, 7:285–297, 1983.
 [42] B. Li and C. A. Fleming. A three dimensional multigrid model for fully nonlinear water waves. Coastal Engineering, 30:235–258, 1997.
 [43] M. S. LonguetHiggins and E. D. Cokelet. The deformation of steep surface waves on water. I. A numerical method of computation. Proc. Roy. Soc. London Ser. A, 350(1660):1–26, 1976.
 [44] H. R. Luth, B. Klopman, and N. Kitou. Projects 13G: Kinematics of waves breaking partially on an offshore bar: LDV measurements for waves with and without a net onshore current. Technical report H1573, Delft Hydraulics, 1994.
 [45] Q. W. Ma, G. X. Wu, and R. Eatock Taylor. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. part 1: methodology and numerical procedure. Int. J. Num. Meth. Fluids, 36(3):265–285, 2001.
 [46] Q. W. Ma, G. X. Wu, and R. Eatock Taylor. Finite element simulations of fully nonlinear interaction between vertical cylinders and steep waves. part 2: numerical results and validation. Int. J. Num. Meth. Fluids, 36(3):287–308, 2001.
 [47] Q. W. Ma and S. Yan. Quasi ALE finite element method for nonlinear water waves. J. Comput. Phys., 212(1):52 – 72, 2006.
 [48] P. A. Madsen and H. A. Schäffer. A review of boussinesqtype equations for gravity waves. In Advances in Coastal and Ocean Engineering, 5:1–95, 1999.
 [49] G. R. Markall, A. Slemmer, D. A. Ham, P. H. J. Kelly, C. D. Cantwell, and S. J. Sherwin. Finite element assembly strategies on multicore and manycore architectures. International Journal for Numerical Methods in Fluids, 71(1):80–97, 2013.
 [50] S. Markidis, J. Gong, M. Schliephake, E. Laure, A. Hart, D. Henty, K. Heisey, and P. Fischer. Openacc acceleration of the nek5000 spectral element code. International Journal of High Performance Computing Applications, 29(3):311–319, 2015.
 [51] H. Mase and J. T. Kirby. Hybrid frequencydomain KdV equation for random wave transformation. In Proc. 23rd International Conference on Coastal Engineering, pages 474–487. ASCE, 1992.
 [52] G. Mengaldo, D. De Grazia, D. Moxey, P. E. Vincent, and S. J. Sherwin. Dealiasing techniques for highorder spectral element methods on regular and irregular grids. J. Comput. Phys., 299:56 – 81, 2015.
 [53] S. B. Nimmala, S. C. Yim, and S. T. Grilli. An efficient 3dfnpf numerical wave tank for virtual largescale wave basin experiments. In In Proc. 31st Intl. Conf. on Ocean, Offshore and Arctic Engineering, 2012.
 [54] A. T. Patera. A Spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput. Phys., 54:468–488, 1984.
 [55] B. T. Paulsen, H. Bredmose, and H. B. Bingham. An efficient domain decomposition strategy for wave loads on surface piercing circular cylinders. Coastal Engineering, 86:57–76, 2014.
 [56] I. Robertson and S. J. Sherwin. Freesurface flow simulation using hp/spectral elements. J. Comput. Phys., 155:26–53, 1999.
 [57] Y.L. Shao and O. M. Faltinsen. A harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics. J. Comput. Phys., 274(0):312 – 332, 2014.
 [58] F. Shi, J. T. Kirby, J. C. Harris, J. D. Geiman, and S. T. Grilli. A highorder adaptive timestepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation. Ocean Modelling, 43–44:36–51, 2012.
 [59] J. Spinneken, V. Heller, S. Kramer, M. Pigott, and A. Viré. Assessment of an advanced finite element tool for the simulation of fullynonlinear gravity water waves. In Proceedings of The Annual International Offshore and Polar Engineering Conference (ISOPE) 2012, 2012.
 [60] V. Sriram, S. A. Sannasiraj, and V. Sundar. Velocity calculation methods in finite element based mel formulation. In Q. W. Ma, editor, Advances in Numerical Simulation of Nonlinear Water Waves, volume 11 of Advances in Coastal and Ocean Engineering, pages 203–244. World Scientific, 2010.
 [61] I. A. Svendsen and I. G. Jonsson. Hydrodynamics of coastal regions. Technical University of Denmark, 2001.
 [62] K. Tanizawa. The state of the art on numerical wave tank. In Proceedings of 4th Osaka Colloquium on Seakeeping Performance of Ships, pages 95–14, 2000.
 [63] M. S. Turnbull, A. G. L. Borthwick, and R. Eatock Taylor. Numerical wave tank based on a transformed finite element inviscid flow solver. Int. J. Num. Meth. Fluids, 42(6):641–663, 2003.
 [64] C.Z. Wang and G.X. Wu. A brief summary of finite element method applications to nonlinear wavestructure interactions. J. Marine. Sci. Appl., 10:127–138, 2011.
 [65] T. Warburton. An explicit construction for interpolation nodes on the simplex. J. Engineering Math., 56(3):247–262, 2006.
 [66] T. Warburton and T. Hagstrom. Taming the CFL number for discontinuous Galerkin methods on structured meshes. SIAM J. Numer. Anal., 46(6):3151–3180, September 2008.
 [67] T. C. Warburton, S. J. Sherwin, and G. E. Karniadakis. Basis functions for triangular and quadrilateral highorder elements. SIAM J. Sci. Comput., 20(5):1671–1695, April 1999.
 [68] G. Wei and J. T. Kirby. Timedependent numerical code for extended boussinesq equations. J. Waterway, Port, Coastal and Ocean Eng., ASCE, 121(5):251–261, 1995.
 [69] J.H. Westhuis. The numerical simulation of nonlinear waves in a hydrodynamic model test basin. PhD thesis, Department of Mathematics, University of Twente, The Netherlands, 2001.
 [70] J.H. Westhuis and A. J. Andonowati. Applying the finite element method in numerically solving the two dimensional freesurface water wave equations. In Proceedings of The 13th International Workshop on Water Waves and Floating Bodies (IWWWFB), 1998.
 [71] G. X. Wu and R. Eatock Taylor. Finite element analysis of twodimensional nonlinear transient water waves. Applied Ocean Res., 16(6):363 – 372, 1994.
 [72] G. X. Wu and R. Eatock Taylor. Time stepping solutions of the twodimensional nonlinear wave radiation problem. Ocean Engineering, 22(8):785 – 798, 1995.
 [73] G. X. Wu and R. Eatock Taylor. The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies. Ocean Engineering, 30(3):387 – 400, 2003.
 [74] V. E. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys., 9:190–194, 1968.
 [75] B. Z. Zhou, G. X. Wu, and B. Teng. Fully nonlinear wave interaction with freely floating nonwallsided structures. Engineering Analysis with Boundary Elements, 50:117 – 132, 2015.
 [76] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. Chapter 6  free surface and buoyancy driven flows. In O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu, editors, The Finite Element Method for Fluid Dynamics (Seventh Edition), pages 195 – 224. ButterworthHeinemann, Oxford, seventh edition edition, 2014.