FNPF-SEM: A parallel spectral element model in Firedrake for fully nonlinear water wave simulations

Abstract

We present a new parallel spectral element solver, FNPF-SEM, for simulating linear and fully nonlinear potential-flow-based water waves and their interactions with offshore structures. The tool is designed as a general-purpose wave model for offshore engineering applications. Built on the open-source Firedrake framework, the new FNPF-SEM model is a computational tool capable of capturing both linear and nonlinear wave phenomena with high accuracy and efficiency, supporting high-order (spectral) finite elements. Additionally, Firedrake provides native support for MPI-based parallelism, enabling efficient multi-CPU distributed computations for large-scale simulations. We demonstrate the capabilities of the high-order model through h- and p-convergence studies, and weak and strong scaling tests. Verification is performed against analytical and numerical solutions and verified against experimental data, including nonlinear high-order harmonic generation and wave interactions with a cylinder and a breakwater. The new FNPF-SEM model offers a numerical framework for simulating wave propagation and wave-structure interactions, with the following key features: (i) the ability to represent complex geometries through flexible, unstructured finite element meshes; (ii) reduced numerical diffusion and dispersion by using high-order polynomial expansions; and (iii) scalability to full- and large-scale simulations over long time periods through a parallel implementation.

Keywords

high-order spectral element method linear and nonlinear potential flow with a free surface water waves wave-structure interactions wave propagation parallel computing

1. Introduction

Simulating the dynamics of water waves remains a significant challenge in computational fluid dynamics (CFD), especially when complex geometries are present in the flow domain and marine or offshore structures are included. However, despite the challenge, accurate modeling of wave propagation and interactions with bathymetry and structures is essential for marine and offshore engineering applications. For example, estimating how the local wave climate and the resulting wave-induced forces impact offshore structures such as large offshore wind turbines Bailey et al. (2015), wave energy devices, and breakwaters. Other cases include metocean studies such as assessments of ocean conditions, Bitner-Gregersen (2015). Many offshore engineering applications must model nonlinear water waves over large regions and over long time scales, requiring effective parallel implementations of these wave models to accurately and efficiently simulate real-world scenarios within acceptable computational times.

Physical accuracy is typically associated with the chosen mathematical model and its representation of the physics. High-fidelity simulations can be based on the Navier-Stokes equations as a model for all fluid flows. However, when solved numerically, these equations come with high computational requirements. For offshore engineering applications, a broad range of free-surface flows can be approximated by fully nonlinear potential flow (FNPF), assuming the fluid is incompressible, inviscid, and irrotational, Dean and Dalrymple (1991). Moreover, if the wave steepness is small, the linear potential flow (LPF) equations are applicable, i.e., a small-amplitude approximation, which reduces the problem’s complexity. One desirable feature of the FNPF formulation is its ability to simulate both offshore and nearshore environments, as wave interactions with the bathymetry are accurately represented, Bihs et al. (2020). When the simulated problem satisfies the assumptions of FNPF theory, a performance improvement of at least two orders of magnitude can be achieved while achieving an accuracy comparable to the full CFD simulations, as shown in Ransley et al. (2019). If the problem adheres further to the assumptions of LPF theory, additional orders of magnitude in computational efficiency can be realized.

Numerically accurate and efficient solutions to FNPF and LPF formulations have been a major research topic for decades; see, e.g., the report by Evans and Newman (1987) from the 1980s and the references herein and the review by Thomas and Dwarakish (2015). Traditionally, the boundary integral method has been popular for wave modeling. To highlight a few of the well-known and recent boundary element-based models (including the traditional low-order boundary element method, the quadratic boundary element method, and the high-order boundary element method), see Liu et al. (2001); Xue et al. (2001); Lee and Newman (2006); Yan and Liu (2011); Harris et al. (2022); Kurnia and Ducrozet (2023); Seixas De Medeiros et al. (2024). Another option are high-order spectral methods, which are very efficient when modeling nonlinear wave propagation and submerged wave-structure interactions, Dommermuth and Yue (1987); Liu et al. (1992); Christiansen et al. (2013); Melander and Engsig-Karup (2024). The FNPF and LPF equations have also been solved using finite difference schemes, Bingham and Zhang (2007); Engsig-Karup et al. (2009); Ducrozet et al. (2010); Amini-Afshar and Bingham (2017, 2018); Bihs et al. (2020), and with harmonic polynomials, Shao and Faltinsen (2014); Tong et al. (2024). Employing the classical second-order accurate finite element method, we have the work by Wu and Eatock Taylor (1994) and many more, e.g., see Le Roux et al. (2000); Wu and Hu (2004); Ma and Yan (2009).

Spectral element methods (SEMs), initially introduced by Patera (1984), combine key advantages of both spectral (polynomial) and finite element methods. SEMs allow for two distinct convergence strategies: (i) h-convergence, where the mesh is refined while maintaining a fixed polynomial order, resulting in algebraic convergence; and (ii) p-convergence, where the mesh is held fixed while increasing the polynomial order of the basis functions, typically leading to spectral (exponential) convergence. The hybrid approach provides efficient, highly accurate numerical solutions by using piecewise polynomial basis functions of arbitrary order with spectral accuracy. Moreover, the SEM can handle complex geometries common in engineering applications because it uses unstructured mesh tessellations. This combination of high accuracy and geometrical flexibility addresses the limitations of traditional numerical approaches, e.g., the difficulty in modeling complex geometry (finite difference) and/or high dispersive and diffusive errors (low-order finite and boundary elements). For further details and applications of the SEM, see the review by Xu et al. (2018) and the textbook by Karniadakis and Sherwin (2005).

Many open-source finite element frameworks that allow for high-order basis functions and parallel computing have been developed over the years. For example, the spectral/hp element framework Nektar++, Cantwell et al. (2015), the arbitrary high-order finite element framework MFEM, Anderson et al. (2021), the deal.II finite element framework, Arndt et al. (2021), or libParanumal, Chalmers et al. (2022). Another computational framework that supports high-order basis functions is Firedrake, Ham et al. (2023). Firedrake is an open-source framework that automates the solution of partial differential equations using finite elements. One of Firedrake’s main features is the unified form language, Alnæs et al. (2014), which allows the symbolic expression of, e.g., weak variational formulations. Moreover, Firedrake runs with Python, uses state-of-the-art PETSc solvers, Balay et al. (1997, 2024), and supports both structured and unstructured meshes, including the most common element types and finite element spaces. Lastly, Firedrake supports fully parallel computation via MPI and can run on large many-core systems. In the recent decade, the Firedrake project has established itself as a robust open-source platform for implementing advanced numerical models.

At DTU, SEM-based FNPF models for simulating free surface water waves have been under active development for more than a decade. This effort began with the work of Engsig-Karup (2006); Engsig-Karup et al. (2006) and later Engsig-Karup et al. (2016), where a solution was proposed to the intrinsic mesh-induced instability associated with finite element discretizations of free surface flows, as previously identified in Robertson and Sherwin (1999). The SEM implementations at DTU are based primarily on the FNPF formulations, Engsig-Karup et al. (2016); Engsig-Karup and Eskilsson (2019); Engsig-Karup et al. (2019); Engsig-Karup and Laskowski (2021); Visbech et al. (2026c), and the LPF approximations, Visbech et al. (2024a,b). This FNPF-SEM approach offers a numerical method that can handle complex geometries and accurately simulate over long time periods using a high-order polynomial basis. The methodology has been extended to encompass the full incompressible Navier-Stokes equations with a free surface in Melander et al. (2024). Achieving low dispersive and diffusive numerical errors arising from the numerical discretization procedure is a key aspect of long-time wave simulations. High-order discretization methods, such as SEM and others, are much more efficient than low-order methods, such as classical p = 1 finite elements, Kreiss and Oliger (1972); Willberg et al. (2012). In addition, Eskilsson et al. (2025) presented a high-performance implementation of an FNPF-based water wave solver using Nektar++. The current capabilities of these FNPF-based SEM models are comparable in functionality, benefiting significantly from the robust, scalable features of modern open-source computational frameworks. Moreover, efficient water wave solutions to the Benney-Luke-type equations, derived using variational principles, are also available via Firedrake, Bokhove and Kalogirou (2016). Furthermore, the Thetis project, an unstructured-grid coastal ocean model, was built using the Firedrake framework, Kärnä et al. (2018). Very recently, steady-state solutions to the incompressible Navier-Stokes equations were presented by Minniti et al. (2025) through a Firedrake implementation.

1.1. Paper contribution

The main aim is to demonstrate the correctness of a new parallel implementation of the aforementioned DTU-based SEM models for solving FNPF and LPF wave problems. The parallel code is implemented in Firedrake, enabling efficient parallel solution of the FNPF equations for large- and full-scale wave problems over long time horizons, ultimately for use in offshore engineering applications. The paper verifies the numerical implementation, referred to as FNPF-SEM, through h- and p-convergence studies and the parallel setup through strong and weak scaling tests. Moreover, the entire model is further verified and validated across various test cases for FNPF and LPF wave propagation and wave-structure interactions.

1.2. Outline of paper

Following the introduction in Section 1, the paper proceeds with the mathematical model in Section 2, which describes the computational domain, its boundaries, and the governing equations. Subsequently, Section 3 outlines the numerical discretization methods, highlighting the chosen spatial and temporal discretizations and key aspects such as wave initialization, generation, and absorption. The results are presented in Section 4, which includes both numerical verification studies (including convergence and scaling analyses) and cases of linear and nonlinear wave propagation and wave–structure interactions. Finally, the paper concludes with a summary of the main findings in Section 5. An appendix containing supplementary material (additional convergence studies, mesh layouts, and simulation results) is included at the back.

2. Mathematical model

Considering a fluid domain, $Ω \in R^{d}$ , where d = 3, and its completely surrounding boundary, Γ = ∂Ω. Here, x = (x, y, z) where x and y span the horizontal plane, and z is the vertical direction defined positively upwards. For the boundary, we have the following decomposition: Γ = ⋃_iΓⁱ for i ∈ {η, b, w, p}. Here, Γ^η is the top free surface boundary defined at z = η in the FNPF setting, where $η = η (x, y, t) : Γ^{η} \times T \mapsto R$ is the free surface elevation defined from z = 0 in the time-domain $T : t \geq 0$ . In the LPF setting, Γ^η is defined at z = 0. Moreover, Γ^b is the bathymetry boundary at the bottom at z = −h, where $h = h (x, y) : Γ^{b} \mapsto R$ is the height of the water column measured around z = 0. At the vertical boundaries, two possibilities exist: Γ^w denotes impermeable walls, and Γ^p denotes periodic boundaries. A conceptual sketch of Ω and its boundaries is shown in Figure 1 for the xz-plane. Remark: Partially or fully submerged structures can be included in this setting, yet this is not the scope of this work.

Figure 1.

Conceptual schematic of the fluid domain, Ω, and its boundaries in the xz-plane.

2.1. Governing equations

Assuming potential flow, i.e., an incompressible, irrotational, and inviscid fluid, the gradient of a scalar velocity potential, $ϕ = ϕ (x, t) : Ω \times T \mapsto R$ , can be used to express the fluid velocities, $u = (u, v, w) = (u (x, t), v (x, t), w (x, t)) : Ω \times T \mapsto R^{d}$ as u = ∇ ϕ where ∇ = (∂_x, ∂_y, ∂_z) is the Cartesian gradient operator in Ω. The spatial mathematical problem is to find ϕ ∈ C²(Ω) by solving the Laplace equation with suitable boundary conditions as

\begin{aligned} \nabla^{2} ϕ & = 0, in Ω, \\ ϕ & = ϕ_{η}, on Γ^{η}, \\ \nabla ϕ \cdot n & = 0, on {Γ^{b}, Γ^{w}}, \\ ϕ |_{Γ^{p}} & = ϕ |_{Γ^{p}}, on Γ^{p}, \end{aligned}

(1)

where

n = (n_{x}, n_{y}, n_{z}) = (n_{x} (x), n_{y} (x), n_{z} (x)) : Γ \mapsto R^{d}

, is the outward-facing normal vector to any boundary of interest. The subscript notation of “η” denotes quantities defined on Γ^η, i.e.,

ϕ_{η} = ϕ_{η} (x, y, t) : Γ^{η} \times T \mapsto R

. This notation becomes useful when we consider the FNPF time governing kinematic and dynamic free-surface conditions in the form of Zakharov (1968). On

Γ^{η} \times T

, we have.

\begin{align} \partial_{t} η & = - \tilde{\nabla} η \cdot \tilde{\nabla} ϕ_{η} + w_{η} (1 + \tilde{\nabla} η \cdot \tilde{\nabla} η), \end{align}

(2)

\begin{align} \partial_{t} ϕ_{η} & = - g η - \frac{1}{2} \tilde{\nabla} ϕ_{η} \cdot \tilde{\nabla} ϕ_{η} + \frac{1}{2} w_{η}^{2} (1 + \tilde{\nabla} η \cdot \tilde{\nabla} η) . \end{align}

(3)

Here $\tilde{\nabla} = (\partial_{x}, \partial_{y})$ is the Cartesian gradient operator only defined on the horizontal free surface, Γ^η, and g = 9.82 [m/s²] is the gravitational acceleration. Given initial conditions for η and ϕ_η at t = 0, this coupled FNPF initial boundary value problem can be solved uniquely for any point in time and space. The classical linearized formulation is obtained by assuming small amplitude-to-length waves, i.e., Γ^η becomes time-invariant and defined at z = 0. Moreover, the two free surface conditions simplify to

\partial_{t} η = w_{η}, and \partial_{t} ϕ_{η} = - g η .

(4)

3. Numerical discretization

We adopt a classical method of lines approach to the governing equations. The spatial and temporal discretizations are discussed below individually. As stated earlier, the mathematical model is implemented in Firedrake, Ham et al. (2023). Firedrake is an open-source computational framework for solving partial differential equations using the finite element method. It provides a high-level interface based on the unified form language, Alnæs et al. (2014), allowing users to express weak variational formulations in a compact, readable form. At the backend, Firedrake automates the generation of optimized C code and supports efficient parallel execution via MPI, making it well-suited for large-scale simulations. Built on PETSc, Balay et al. (2024), it includes state-of-the-art solvers and supports a wide range of element types and mesh structures, including high-order spectral elements on both structured and unstructured meshes.

3.1. Spatial discretization

We perform a discrete tesselation of the spatial domain, Ω, and boundary, Γ^η, using non-overlapping finite elements as

D ≃ T_{D} = ⋃_{n = 1}^{N_{D}} D_{n}, for D \in {Ω, Γ^{η}},

where

N_{D}

is an arbitrary number of elements, such that

D_{n}

is the n’th element in the tesselation

T_{D}

. The element types and meshes are highlighted in the following section. On the two tesselations, we introduce

C^{0} (D)

function spaces of piecewise polynomial functions as:

V_{D} = {ν_{D} \in C^{0} (D); \forall n \in {1, \dots, N_{D}}, ν_{D} |_{D_{n}} \in P^{p}},

where

P^{p}

denotes the polynomial space of nodal Lagrange polynomials of degree at most p on elements with Gauss-Lobatto-Legendre (GLL) node distribution, Rønquist and Patera (1987).

3.1.1. Weak variational formulations

Next, we state the weak forms of (1) and (2) through (3), respectively. The linearized forms in (4) are omitted for brevity. The linear variational problem for the Laplace problem is to find $ϕ \in V_{Ω}$ such that $\forall ν_{Ω} \in V_{Ω}$ , with

a [ϕ; ν_{Ω}] = b [ν_{Ω}],

(5)

where the bi-linear and linear forms are derived by multiplying by a test function, ν_Ω, integrating over Ω, and applying integration by parts as

\begin{array}{l} a [ϕ; ν_{Ω}] & = {⟨\nabla ϕ, \nabla ν_{Ω}⟩}_{Ω}, \\ b [ν_{Ω}] & = 0, \end{array}

with ⟨a,b⟩_S = ∫_S a ⋅ b dS denoting the inner product over the set S.

For the two free surface conditions, we define a nonlinear variational problem to find ${\partial_{t} η, \partial_{t} ϕ_{η}} \in V_{Γ^{η}}$ such that $\forall ν_{Γ^{η}} \in V_{Γ^{η}}$ , we solve

F_{k} [\partial_{t} η; ν_{Γ^{η}}] = 0, and F_{d} [\partial_{t} ϕ_{η}; ν_{Γ^{η}}] = 0,

(6)

where the nonlinear variational forms are derived by multiplying by a test function,

ν_{Γ^{η}}

, integrating over Γ^η, and moving everything to the left-hand side. Ultimately, this gives

\begin{aligned} F_{k} [\partial_{t} η; ν_{Γ^{η}}] = {⟨\partial_{t} η, ν_{Γ^{η}}⟩}_{Γ^{η}} + \\ {⟨\tilde{\nabla} η \cdot \tilde{\nabla} ϕ_{η} - w_{η} (1 + \tilde{\nabla} η \cdot \tilde{\nabla} η), ν_{Γ^{η}}⟩}_{Γ^{η}}, \end{aligned}

and

\begin{aligned} F_{d} [\partial_{t} ϕ_{η}; ν_{Γ^{η}}] = {⟨\partial_{t} ϕ_{η}, ν_{Γ^{η}}⟩}_{Γ^{η}} + \\ {⟨g η - \frac{1}{2} \tilde{\nabla} ϕ_{η} \cdot \tilde{\nabla} ϕ_{η} - \frac{1}{2} w_{η}^{2} (1 + \tilde{\nabla} η \cdot \tilde{\nabla} η), ν_{Γ^{η}}⟩}_{Γ^{η}} . \end{aligned}

3.1.2. Initial meshes, extrusion, and mesh updating

Firedrake has native support for extruded meshes, McRae et al. (2016); Bercea et al. (2016). The mesh construction for this model proceeds in three steps: First, a triangular mesh for the free surface, $T_{Γ^{η}}$ , is generated. This is done using Firedrake’s built-in mesh utilities for simple rectangular domains or via Gmsh for more complex geometries, Geuzaine and Remacle (2009). Second, the surface mesh is vertically extruded in the z-direction into N_z layers to form the three-dimensional fluid domain mesh, $T_{Ω}$ . These two meshes, referred to as the ”initial meshes” $(T_{Γ^{η}}^{i}$ and $T_{Ω}^{i})$ , consist entirely of prismatic elements, and can be unstructured to adapt the geometric features in the horizontal (x, y)-plane and structured in the z-direction. This means that $T_{Ω}$ has quadrilateral vertical boundary meshing and triangular top and bottom boundary meshing. Third, in the FNPF setting, a vertical translation is applied to each column of points to track the evolving free-surface curvature. This ensures an accurate geometric representation of the time-dependent fluid domain, Ω, defined by the free-surface elevation, η. Consequently, mesh updates are performed at each discrete time stage or step to reflect the dynamic boundary. In Figure 2 (lower), the updated curvilinear mesh can be seen. Note that elements in Firedrake (potentially of high order and curved) are supported through iso-parametric mappings, Rathgeber et al. (2016). Remark: A similar approach is used to modify the mesh due to non-flat bathymetries. Currently, submerged and floating structures are not directly supported in the model due to the mesh extrusion; however, work is underway to address this.

Figure 2.

Point distribution (indicated by different colors) for a finite FNPF wave tank distributed over N_c = 5 cores on a mesh with (N_x, N_y, N_z) = (20, 4, 4) elements of order p = 5.

We want to move information and data (when running in parallel) between the free surface mesh, $T_{Γ^{η}}$ , and the corresponding top free surface points of the fluid domain mesh, $T_{Ω}$ . Using extruded meshes in Firedrake, it is guaranteed that the parallelized partition of the two meshes coincides (in the x- and y-coordinates) as seen in Figure 2. This feature makes it possible to compute a mapping, $m : T_{Ω} \mapsto T_{Γ^{η}}$ , that extracts the free surface values, $f_{Γ^{η}}$ , from any fluid domain variable, f _Ω, as $f_{Γ^{η}} = f_{Ω} [m]$ , i.e., a simple index/mask operation between two vectors.

3.2. Temporal discretization and modal filtering

Time integration is performed using a standard explicit 4th-order, 4-stage, Runge-Kutta scheme (ERK4). The time step size, Δt, is kept constant throughout the simulation and is computed based on a Courant-Friedrichs-Lewy (CFL) condition: $Δ t = C_{CFL} Δ x_{\min} u_{\max}^{- 1}$ . Here, C_CFL = 0.95, Δx_min denotes the minimum spacing between two grid points on Γ^η, and u_max represents the maximum velocity encountered in the simulation. Remark: By choosing the ERK4, the temporal discretization becomes high-order accurate, and this way complements the high-order spatial discretization. It is noted that the time step size is mostly physics (depth)-restricted rather than restricted by the horizontal spatial resolution, as is common in typical hyperbolic PDE solvers, as discussed in (Engsig-Karup et al., 2016: Section 4.1). Hence, the scheme may benefit from incorporating adaptive time-step control for improving the overall efficiency of the simulations.

Some form of stabilization is typically required to accurately model strongly nonlinear wave phenomena. In Engsig-Karup et al. (2016), the instability issue initially highlighted by Robertson and Sherwin (1999) was addressed using structured elements aligned in the vertical direction. This was coupled with an over-integration (collocation) approach for the nonlinear free surface terms, Kirby and Karniadakis (2003), and a mild modal (spectral) filter applied to the free surface quantities, Hesthaven and Kirby (2008). In Firedrake, the integration order of individual integral terms can be explicitly specified, facilitating the over-integration strategy suitable for reducing aliasing instabilities of nonlinear wave propagation, Engsig-Karup et al. (2016). A modal filter was implemented within this framework by developing a custom PyOP2 kernel written in C, executed via the par_loop interface. This allows the filter to be applied in parallel across all free-surface elements using the interior point values. Specifically, the kernel operates element-wise, applying the local modal filter as defined in Engsig-Karup et al. (2016): Suppose g_n(x) is a function defined on a free surface element $Γ_{n}^{η}$ . On each element, we have n_p = (p + 1) (p + 2)/2 nodal points denoted by x _n = (x_n, y_n). We map the solution to the modal space as

{\hat{g}}_{n} (f) = V^{- 1} g_{n} (x_{n}),

where

V \in R^{n_{p} \times n_{p}}

is a conforming Vandermonde matrix defined as

V_{i j} = {\hat{ψ}}_{j} (r_{i})

, where r _i = (r_i, s_i), is coordinates in a reference domain,

R = {r = (r, s) | (r, s) \geq - 1; r + s \leq 0}

and

\hat{ψ} (r)

is a modal basis function. See Engsig-Karup et al. (2016) for more details. Using a filter matrix,

F \in R^{n_{p} \times n_{p}}

, filtering can be applied to the highest modes

{\hat{g}}_{n, f i l t e r e d} (f) = F {\hat{g}}_{n} (f) .

Then, the filtered solution can be mapped back to nodal space

g_{n, f i l t e r e d} (x_{n}) = V {\hat{g}}_{n, f i l t e r e d} (f) .

This filtering approach can be applied to any free surface quantity, e.g., η, ϕ_η, or w_η, as needed for nonlinear waves.

3.3. Wave initialization, generation, and absorption

For periodic wave tank simulations, the model is initialized with η and ϕ_η, either from the linear wave theory in the LPF setting or from an accurate stream function solution in the FNPF setup, Fenton (1988). Linear (small-amplitude) waves can be fully characterized by the non-dimensional wave number, kh, given a specified wavelength, L, or water depth, h. In contrast, nonlinear (finite-amplitude) stream function waves also require the specification of a relative wave steepness, ɛ/ɛ_max, to quantify their nonlinearity. Here, ɛ = H/L denotes the wave steepness based on the wave height H, and ɛ_max is the maximum steepness attainable before wave breaking for a given wave, Fenton (1990). The initial conditions are set to zero for simulations in numerical wave tanks of finite size. Wave generation and absorption are handled using classical relaxation zones, Larsen and Dancy (1983), of 2-3 wavelengths in size. The wave fields, η and ϕ_η, are then modified at each time step as a post-processing operation as

f = (1 - C_{r}) f + C_{r} f_{true},

where f ∈ {η, ϕ_η} and

C_{r} = C_{r} (x, y) : Γ^{η} \mapsto R

is the relaxation function, taking the value of zero in the target zone and smoothly goes to unity in the generation and absorption zone. In the generation zone, f_true is the true actual solution and zero outside this zone. To avoid impulsive behavior at the onset of the simulation, the generation zone is gradually activated from 0% to 100% over five wave periods.

3.4. Entire computational workflow

To support the understanding of how the entire computational workflow is outlined, a step-by-step procedure is shown below:

Step 1 – Setup: Build initial meshes, $T_{Γ^{η}}^{i}$ and $T_{Ω}^{i}$ , and compute pre-computable quantities, e.g., the map, m.

Step 2 – Initial conditions: Impose initial conditions on η and ϕ_η at t = t_start.

Step 3 – Time integration: Iterate Step 4 (below) over time steps/stages until t = t_end.

Step 4.1 – Mesh-update: Update $T_{Ω}$ according to η by shifting the vertical coordinate z and using the map, m, to transfer information.

Step 4.2 – Boundary conditions: Impose direct Dirichlet boundary conditions using ϕ_η.

Step 4.3 – Laplace problem: Solve the Laplace problem for the velocity potential, ϕ.

Step 4.4 – Gradient recovery: Compute the vertical velocity, w, from ϕ via gradient recovery.

Step 4.5 – Vertical free surface velocity: Extract the vertical free surface velocity, w_η, from w via the map, m.

Step 4.6 – Free surface conditions: Evaluate the two free surface conditions using η, ϕ_η, and w_η.

Step 4.7 – Relaxation and filtering: Modify wave solution depending on generation and absorption zones and apply filtering as needed.

Step 5 - Post-processing: Plot results, etc.

4. Results

In the following, we present numerical results for the model’s verification and validation. All simulations are performed on the HPC system at the DTU Computing Center, Technical University of Denmark, using the hardware configurations listed in Table 1.

Table 1.

Hardware configurations.

Parameter(s)	Value
Total nodes, N_n, available	15
Total cores, N_c, available	300
Cores per node, N_c/N_n	20
Memory per node	125 GB
CPU	Intel Xeon E5-2660 v3
CPU frequency	2.60 GHz

As for the numerical solver strategies, we solve the Laplace problem in (1) following Rathgeber et al. (2016) by using a CG solver with the HYPRE BoomerAMG algebraic multigrid preconditioner, Hestenes and Stiefel (1952); Falgout et al. (2006), using a relative and absolute tolerance of 10⁻⁶ and 10⁻¹⁵, respectively. For the evaluation of the two nonlinear free surface conditions in (2) and (3), we employ the default solver using a GMRES strategy with an ILU as preconditioner and a relative and absolute tolerance of 10⁻⁵ and 10⁻¹⁵, respectively. For this study, this leads to an efficient solver strategy for the cases considered, cf. the analysis of iterative solver efficiency for FNPF modeling reported in Engsig-Karup (2014).

4.1. Model verification

During model verification, various timings will be presented. Those timings are an average for the number of cores, N_c, used in the test. When using N_c ≤ 20, a complete node is allocated to the run.

4.1.1. Convergence h/p-studies

Investigating and confirming the expected convergence rates of the proposed model constitutes the first step in the verification process. For this, we perform h- and p-convergence studies; that is, mesh refinement under a fixed polynomial order of the basis functions, and an increase in the polynomial order on fixed meshes, respectively. For the analysis, a periodic domain is considered, featuring nonlinear waves for all combinations of kh = {1, 3, 6} and ɛ/ɛ_max = {10%, 50%, 90%}. The error is evaluated in the ∞-norm of the vertical free surface velocity, $| w_{η}^{true} - w_{η} |_{\infty}$ , where $w_{η}^{true}$ denotes the reference (true) solution and w_η is the corresponding numerical solution. All wave solutions are defined over a domain of fixed length L = 1 [m], implying L_x = L_y = L. The depth is varied depending on the specific value of kh.

The results of the h-refinement convergence study are presented in Figure 3, using basis functions of order p = {1, 2, 3, 4, 5}. The parameter h_max denotes the maximum length of the element in the computational mesh. Due to the gradient recovery of w from ϕ, it is well known that an order of algebraic convergence is typically lost, Engsig-Karup et al. (2016), leading to an expected convergence rate of $O (h_{\max}^{p})$ . This expected behavior is clearly observed in the results. Furthermore, increasing the wave’s nonlinearity is shown to increase the problem’s complexity, requiring finer meshes to attain comparable levels of accuracy.

Figure 3.

Convergence study for h-refinement for nonlinear waves in all combinations between kh = {1, 3, 6} and ɛ/ɛ_max = {10%, 50%, 90%}.

The p-convergence study is presented in Appendix C, in Figure 19, where spectral, i.e., exponential, convergence is observed. The test is conducted on six different meshes with increasing resolution. Specifically, Mesh 1 has h_max = 0.25, Mesh 2 has h_max = 0.125, Mesh 3 has h_max = 0.0625, etc. As observed previously, increasing nonlinearity makes the problem more challenging, requiring higher resolution to achieve comparable accuracy. Remark: Similar analyses have been performed for the LPF model, which exhibits improved convergence behavior. This is an expected outcome for a linearized model, and therefore, the results are omitted here for the sake of brevity.

4.1.2. Performance profiling

It is well established that solving the Laplace problem constitutes the primary computational bottleneck in both LPF and FNPF wave models. Distributing the finite element domain across multiple cores inevitably affects performance. To this end, we profile the three most computationally intensive routines involved in time integration: (1) LaplaceSolve, which includes imposing boundary conditions, solving the linear system of equations, and performing gradient recovery (Step 4.2, 4.3, 4.4, 4.5); (2) EvaluateRHS, responsible for evaluating the two free surface boundary conditions (Step 4.6); and (3) LaplaceUpdate, which updates the fluid domain and associated functions (Step 4.1).

Profiling is carried out on N_c = {1, …, 300} cores using two different meshes of polynomial order p = 4 with fixed vertical resolution N_z = 4, as detailed in Table 2. Mesh 2 is four times larger than Mesh 1, both in terms of the number of elements and degrees of freedom (DoF).

Table 2.

Performance profiling meshes with fixed N_z = 4 and p = 4.

Mesh	N _Ω	$N_{Γ^{η}}$	$DoF (T_{Ω})$	$DoF (T_{Γ^{η}})$
1	32,768	8,192	1,114,112	65,536
2	131,072	32,768	4,456,448	262,144

Figure 4 presents the core-averaged timing for each time stage across the three routines. A substantial reduction in time per stage is observed; the increase is from N_c = 1 to approximately N_c = 20. Beyond this point, only the LaplaceSolve routine continues to benefit from further parallelization, while EvaluateRHS and LaplaceUpdate begin exhibiting constant runtimes. Moving to a finer mesh (Mesh 2) increases the computational time for both LaplaceSolve and EvaluateRHS, as expected. In contrast, LaplaceUpdate remains largely unaffected by mesh resolution, suggesting that its performance is dominated by communication overhead. The stagnation in EvaluateRHS scalability is attributed to the relatively small size of the free surface mesh, which limits the parallel efficiency of this part. However, the overall solver remains scalable as long as LaplaceSolve is the dominant cost factor. Moreover, comparing the results of Mesh 1 and 2 shows that increasing the problem size enables scalable behavior across more cores, as LaplaceSolve is the dominant cost factor for longer runs.

Figure 4.

Timings for routines in the FNPF model for Mesh 1 (blue) and Mesh 2 (red), using a polynomial order p = 4.

Figure 5 shows the relative percentage distribution of the total run time across the three routines. At low core counts (e.g., N_c ≤ 5), LaplaceSolve dominates the computational cost. However, as N_c increases, its relative share decreases because the remaining routines incur a constant cost. For larger problems (as in Mesh 2), this imbalance is somewhat mitigated. Potential optimizations include reducing the core count assigned to the free-surface evaluation, as done in Eskilsson et al. (2025), or implementing more efficient iterative solvers for the nonlinear free-surface problem. These optimizations are left for future work. Remark: The LPF formulation does not exhibit these performance limitations, as it avoids both the domain update step to solve (5) and the nonlinear variational problem in (6), which reduces to a simple function evaluation in the linearized setting.

Figure 5.

Percentage distribution of computational cost across routines for Mesh 1 (blue) and Mesh 2 (red), using p = 4.

Based on this code performance profiling, it is confirmed that the Laplace solver is the primary computational bottleneck. Consequently, the following strong- and weak-scaling analyses focus exclusively on this part of the code. Moreover, the convergent behavior showcased in the previous section has been confirmed on multicore hardware setups as well; however, it is omitted for conciseness.

4.1.3. Parallel performance: Strong scaling

To evaluate how runtime performance is affected by increasing the number of cores for a fixed problem size (i.e., a constant total workload), we perform strong-scaling tests of the FNPF model. These tests assess both parallel efficiency and speedup as a function of the number of compute nodes, N_n, where each node comprises 20 cores. As a metric for strong scaling efficiency, we define γ_s as

γ_{s} = \frac{T^{b} N_{n}^{b}}{T_{N_{n}} N_{n}},

(7)

where T^b and

N_{n}^{b}

denote the runtime and node count of the baseline simulation, respectively, and

T_{N_{n}}

is the runtime using N_n nodes. In this study, the baseline corresponds to a single node configuration

(N_{n}^{b} = 1)

. Strong scaling experiments are conducted using the three mesh configurations listed in Table 3, with polynomial orders p = {3, 4, 5} and a fixed vertical resolution of N_z = 4. The tests are executed on up to N_n = 15 nodes, equivalent to N_c = 300 cores, using a nonlinear periodic domain with kh = 2π.

Table 3.

Strong scaling meshes with N_z = 4.

Mesh	N _Ω	$N_{Γ^{η}}$	p	$DoF (T_{Ω})$	$DoF (T_{Γ^{η}})$
1	32,768	8,192	3	479,232	36,864
			4	1,114,112	65,536
			5	2,150,400	102,400
2	131,072	32,768	3	1,916,928	147,456
			4	4,456,448	262,144
			5	8,601,600	409,600
3	262,088	65,522	3	3,833,037	294,849
			4	8,910,992	524,176
			5	17,199,525	819,025

Figure 6 shows the strong scaling efficiency using (7) for the three meshes. As expected, both increased mesh refinement and higher polynomial order improve the model’s parallel efficiency. This reflects better utilization of computational resources for larger problem sizes. The corresponding speedup curves are presented in Figure 7, alongside an ideal speedup line. The results indicate that the model approaches optimal speedup behavior as the problem size increases. For smaller problems, however, speedup plateaus early, revealing limitations in parallel efficiency due to insufficient workload per core.

Figure 6.

Strong scaling efficiency test. Mesh 1 (blue), Mesh 2 (red), and Mesh 3 (orange) with polynomial orders p = {3, 4, 5}.

Figure 7.

Speedup (compared to N_n = 1 as baseline) for the strong scaling test. Mesh 1 (blue), Mesh 2 (red), and Mesh 3 (orange) with polynomial orders p = {3, 4, 5}.

4.1.4. Parallel performance: Weak scaling

To assess how the model scales with a constant workload per core, we conducted a weak-scaling analysis. The test domain is periodic, featuring nonlinear waves characterized by kh = 2π and L = 1 [m]. The resolution is fixed in both the y- and z-directions with N_y = N_z = 4. The simulations are run on N_n = {1, …, 15} nodes corresponding to N_c = {20, …, 300} cores. To ensure a constant workload per core, the domain length is set as L_x = N_cL, which means that each core is responsible for simulating one wave within a domain of length L. The mesh configurations used for the tests are detailed in Table 4, all using polynomial order p = 5. The largest simulation, executed on 300 cores, corresponds to a problem size of

DoF (T_{Ω}) = 80,640,000

. The weak scaling efficiency is defined as

γ_{w} = \frac{T^{b}}{T_{N_{n}}} .

(8)

Table 4.

Weak scaling meshes with N_y = N_z = 4 and polynomial order p = 5. All reported quantities are per core, N_c.

Mesh	N _Ω	N _x	$DoF (T_{Ω})$	$DoF (T_{Γ^{η}})$
1	128	4	8,400	400
2	256	8	16,800	800
3	512	16	33,600	1,600
4	1,024	32	67,200	3,200
5	2,048	64	134,400	6,400
6	4,096	128	268,800	12,800

Figure 8 presents the run times for the various mesh configurations. As expected, increasing resolution increases run time. However, more importantly, the runtime remains nearly constant as the number of cores increases, confirming the model’s weak scaling behavior. This trend is further analyzed in Figure 9, which shows the weak scaling efficiency using (8). The results indicate that higher mesh resolutions lead to improved parallel efficiency, as the computational workload per core becomes more substantial relative to communication overhead. Remark: Additional tests have been conducted for other polynomial orders, showing similar trends, as well as weak scaling tests for fixed meshes and increasing polynomial order. These results have been omitted for brevity.

Figure 8.

Weak scaling timings on six meshes with p = 5.

Figure 9.

Weak scaling efficiency on six meshes with p = 5.

4.2. Simulations

For all the following cases, the simulations are run on a full node (N_c = 20 cores).

4.2.1. Generation of high-order wave harmonics over a submerged bar

First, we consider the classical bar test, which can be compared with the physical experimental laboratory data by Beji and Battjes (1994), denoted as BB94 in the following. Note that this is a wave flume test, thus two-dimensional; however, we extend the domain in the y-direction to match our three-dimensional model. Mildly nonlinear waves are generated and propagated over a submerged bar, which generates high-order wave harmonics due to nonlinear shoaling effects (steepening and shortening of the waves). The domain is sketched out conceptually in Figure 1, where the specifications of the domain are: Generation zone (0 ≤ x ≤ 8), target zone with the submerged bar (8 ≤ x ≤ 30), and absorption zone (30 ≤ x ≤ 38). The width of the domain is 0 ≤ y ≤ 1. All units are [m]. The depths before and after the bar are h = 0.4 [m] and h = 0.1 [m] above the bar. The incline before and decline after the bar are 1:20 and 1:10, respectively. The input wave characteristics are set as kh = 0.6725 [-] and ɛ/ɛ_max ≈ 7% given a wave period of T = 2.018 [s]. For the simulation, we used a free-surface mesh of 7600 elements, extruded with N_z = 4 elements in the vertical direction, yielding 22,800 prismatic elements to represent Ω. Ultimately, this yields problem sizes of 810,693 DoF and 62,361 DoF for Ω and Γ^η, respectively, when using basis functions of order p = 4. The simulation is run for 25 wave periods.

A contour plot of the free surface elevation, η, at the final time is shown in Figure 10. Here, the high-order harmonic is observed to occur when the wave reaches the submerged bar. In addition, a harmonic decomposition of the wave signal is carried out and presented in Figure 11. The harmonic analysis is performed by fitting the time series, both experimental and numerical, using a least-squares approach to a sum of sinusoidal functions with argument 2πf_nt, where f_n = n/T and n = {1, …, 4}. The figure illustrates how waves with dominant first-order components, along with minor second-order contributions, propagate toward the bar. Over and beyond the bar, second-order components begin to dominate the wave dynamics, with contributions from the first, third, and fourth harmonics. Finally, in the absorption zone, all wave energy is dissipated. The comparison between the numerical solution (FNPF) and the experimental data (BB94) shows good visual agreement.

Figure 10.

Contour of the free surface elevation at t = 25T.

Figure 11.

Harmonic analysis of the numerical simulation (FNPF) and experimental data (BB94) for the most dominant harmonic modes in the bar test. A₁ = H/2 = 0.01 [m].

4.2.2. Linear and nonlinear interactions with a vertical cylinder

Next, we investigate wave diffraction and the maximum free surface elevation, i.e., wave run-up, on a vertical surface-piercing cylinder subjected to both linear and nonlinear waves. The computational domain is symmetric about the y = 0 plane, and the outer bounds are defined by 0 ≤ x ≤ 12L and 0 ≤ y ≤ 8L. Wave generation and absorption zones are located in the regions 0 ≤ x ≤ 2L and 10L ≤ x ≤ 12L, respectively. The vertical cylinder, of radius R, is centered at (x, y) = (6L, 0). Each simulation is run for 25 wave periods. The horizontal mesh is fully unstructured and graded to provide higher resolution near the cylinder wall boundary Γ^cylinder, with approximately eight elements per radius R. In contrast, far-field boundaries are meshed more coarsely, using roughly two elements per wavelength L. The vertical discretization uses N_z = 4 layers. The cylinder boundary is currently approximated using affine geometry, with a higher-order curvilinear representation left for future development. All relevant parameters for the domain, mesh, and wave conditions are summarized in Table 5 for both the linear and nonlinear cases. In Appendix A, Figure 15, the unstructured mesh used to obtain the linear solution can be seen.

Table 5.

Parameters for cylinder simulations.

	Linear	Nonlinear
Test case (TC)	1	2
R	0.5 (m)	0.1625 (m)
L	1 (m)	2.73 (m)
h	1/(2π) (m)	0.45 (m)
kh	1 (-)	1.036 (-)
H	0.01 (m)	0.053 (m)
ɛ/ɛ_max	0%	19%
p	4	5
N _Ω	8,824	15,048
$N_{Γ^{η}}$	2,206	3,762
$DoF (T_{Ω})$	305,133	1,000,251
$DoF (T_{Γ^{η}})$	17,949	47,631

The results from both simulations are shown in Figure 12, presented as the non-dimensional maximum free-surface elevation at the cylinder boundary, η_m/(2H). Here, η_m is

η_{m} (x_{i}, y_{i}) = \max_{t} η (x_{i}, y_{i}, t), (x_{i}, y_{i}) \in Γ^{cylinder} .

Figure 12.

Maximum non-dimensional free surface elevation on the cylinder, η_m/(2H), for the two simulations comparing numerical, experimental, and analytical solutions.

For the linear test case (TC 1), the numerical results, denoted LPF-1, are compared with the analytical predictions of linear diffraction theory by MacCamy and Fuchs (1954), denoted MF54-1. For the nonlinear test case (TC 2), the numerical results of the FNPF simulation, denoted FNPF-2, are compared with the digitized experimental data of Kriebel (1992), denoted K92-2, which were digitized for comparison. We also include equivalent linear results for this configuration, denoted MF54-2.

Figure 12 shows good agreement between the numerical and analytical results for the linear case (LPF-1 vs MF54-1), particularly in the region downstream of the cylinder. For the nonlinear test case, the numerical solution (FNPF-2) aligns well with the experimental data (K92-2), consistent with previous results reported in the literature (e.g., see Ducrozet et al. (2010)). Notably, nonlinear effects (compared to MF54-2) manifest as increased run-up at 0° and 180°, and a reduced run-up at intermediate angles, around 45° and 315°.

4.2.3. Wave interactions with a V-shaped breakwater

As a final verification test case, we consider wave interactions with a breakwater structure, following the setup described in Chang et al. (2012), for incident waves approaching at angles of θ = {0°, 30°, 60°}. The computational domain is defined as follows: From the origin, (x, y) = (0, 0), two breakwater arms extend a length of b = 76.2 [m] and form an internal angle of 2β = π/3 [rad] (60°). We refer to the upper and lower arms as the upper breakwater arm (UBA) and lower breakwater arm (LBA). At the ends of the arms, inward-facing curved end-walls are placed, each spanning an angle of 2γ = 2π/9 [rad] (40°). While the original study by Chang et al. (2012) treated the breakwaters as infinitely thin, we adopt a more physically realistic representation by assigning them a finite thickness of l = 1 [m].

The computational domain spans the region −L₂ − 2b ≤ (x, y) ≤ 2b + L₂, completely enclosing the target region defined as −2b ≤ x ≤ 2b and −L₂ − 2b ≤ y ≤ 2b + L₂. Here, L₂ = 260 [m] $\approx 2 L$ , with L denoting the wavelength of the incoming wave. A generation zone is defined for −L₂ − 2b ≤ x ≤ −2b, and an absorption zone for 2b ≤ x ≤ 2b + L₂. The water depth is set to h = 13.4 [m] throughout. For the spatial discretization, a fully unstructured horizontal mesh is employed. The mesh is coarsened toward the domain ends (approximately one element per wavelength, L) and refined near the breakwaters (approximately 10 elements per breakwater arm length, b). This results in 754 triangular elements on the free surface, $T_{Γ^{η}}$ , which are extruded into four vertical layers, yielding a total of 3,016 prismatic elements in the fluid domain, $T_{Ω}$ . Using a polynomial order of p = 5, the resulting discrete system contains 202,020 degrees of freedom (DoF). See Figure 13 for $T_{Γ^{η}}$ .

Figure 13.

Free surface mesh for θ = 0° near the breakwater.

The wave input are from Chang et al. (2012), with a non-dimensional wavenumber kh = 0.65 [-], corresponding to L = 129.5 [m]. The linear wave height is set to H = 1.295 [m], and the simulations are run for 10 periods. For each incident angle, θ = {0°, 30°, 60°}, we compute the variation in η at each grid point (x_i, y_i) as

η_{v} = η_{v} (x_{i}, y_{i}) = \max_{t} η (x_{i}, y_{i}, t) - \min_{t} η (x_{i}, y_{i}, t) .

The corresponding non-dimensional quantity η_v/H is visualized as contour plots in Appendix B (Figures 16 –18). In Figure 14, we present η_v/H along the UBA and LBA for the three wave angles, comparing our numerical results (LPF-θ) with digitized reference data from Chang et al. (2012) (C12-θ). Despite small differences arising from the use of finite-thickness structures, domain discretization, and minor wave reflections at the boundaries, the numerical solutions show good qualitative agreement with the reference results. This highlights that the incident angle affects the wave elevation in the vicinity of the breakwater structure, e.g., which may be of interest for sea-state or load estimations on the structure. Added to the figure, we also show the results from the FNPF model (FNPF-θ) with the same wave input parameters; however, with an incident wave of ɛ/ɛ_max = 18%. From the figure, the added nonlinearity is visually evident in the maximum variance of the free surface elevation on the structure, especially on the LBA side. For θ = 60° on the UBA side, only a minor change can be observed; however, that is also the side completely sheltered from the incident waves.

Figure 14.

Breakwater solution (non-dimensional free surface variation) on the UBA and LBA for θ = {0°, 30°, 60°} using proposed linear (LPF-θ) and nonlinear (FNPF-θ) numerical model and results from Chang et al. (2012) (C12-θ).

5. Conclusion

In this paper, we presented a new parallel implementation of a linear and fully nonlinear potential flow model, FNPF-SEM, to simulate wave propagation and wave-structure interaction. The model employs high-order (spectral) finite elements within the open-source Firedrake framework, leveraging its native MPI-based parallelism to handle large-scale computations efficiently. We detailed the mathematical formulation, including the governing equations, domain definitions, and numerical discretization. Parallel performance was assessed through code profiling for strong- and weak-scaling studies, demonstrating good parallel efficiency, particularly for large-scale simulations. Profiling confirmed that the Laplace problem constitutes the dominant computational cost. The numerical accuracy of the implementation was verified through convergence studies, confirming the expected algebraic and spectral convergence under mesh refinement and increased polynomial order, respectively. Finally, the model was compared against analytical solutions and experimental data for a range of wave-propagation and wave-structure-interaction scenarios, including high-order harmonic generation over a submerged bar and linear and nonlinear wave interactions with a vertical cylinder and a V-shaped breakwater. These results confirm the accuracy and applicability of the proposed framework. In ongoing work, the tool will be expanded to handle regional-scale wave propagation for accurate sea-state estimation and to consider areas with irregular coastlines, exploiting the advantages of unstructured meshes.

Footnotes

Acknowledgements

The research was carried out at the Technical University of Denmark (DTU) in the Department of Applied Mathematics and Computer Science, with computational resources provided by the DTU Computing Center (DCC). Lastly, the authors express their gratitude to the Firedrake team led by Professor David Ham for their support.

ORCID iDs

Jens Visbech

Anders Melander

Allan Peter Engsig-Karup

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work partly contributes to the activities of the research project: “A new digital twin concept for floating offshore structures” supported by COWIfonden (Grant no. A-165.19) and partly the PhD-project of JV: “New Advanced Simulation Techniques for Wave Energy Converts”.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data Availability Statement

The computational model developed in this study is publicly available on GitHub. The repository contains the full implementation, including source code, documentation, and example scripts necessary to reproduce the results presented in this paper. Additional information or clarification can be obtained from the corresponding author upon reasonable request. Link to GitHub: *

Author biographies

Jens Visbech is a PhD student in Scientific Computing at the Department of Applied Mathematics and Computer Science, Technical University of Denmark. His research focuses on the development of high-order finite element methods for wave–structure interaction problems, with applications to renewable energy systems such as floating offshore wind turbines and wave energy converters.

Anders Melander is a postdoc from the Section for Scientific Computing, Department of Applied Mathematics and Computer Science, Technical University of Denmark. He received his PhD from the same affiliation. His research focuses on high-order numerical methods for nonlinear free-surface flows.

Allan Peter Engsig-Karup is an Associate Professor of Computational Mathematics in the Section for Scientific Computing, Department of Applied Mathematics and Computer Science, Technical University of Denmark. He is affiliated with the Danish Center for Applied Mathematics and Mechanics (DCAMM) and the Center for Energy Resources Engineering (CERE). His research interests relates to modern scientific computing and include high-performance computing, multi-device GPU acceleration, numerical analysis, and the development of efficient and scalable algorithms. His work focuses on researching advanced numerical methods and developing new state-of-the-art simulation tools based on high-order numerical methods, with applications in computational hydrodynamics (free surface flows), computational acoustics, computational fluid dynamics, reservoir simulation and physics-informed machine learning with applications related to energy resource engineering.

Appendix

References

Alnæs

Logg

Ølgaard

, et al. (2014) Unified form language: a domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software 40(2): 1–37. https://doi.org/10.1145/2566630

Amini-Afshar

Bingham

(2017) Solving the linearized forward-speed radiation problem using a high-order finite difference method on overlapping grids. Applied Ocean Research 69: 220–244. https://doi.org/10.1016/j.apor.2017.11.001

Amini-Afshar

Bingham

(2018) Pseudo-impulsive solutions of the forward-speed diffraction problem using a high-order finite-difference method. Applied Ocean Research 80: 197–219. https://doi.org/10.1016/j.apor.2018.08.017

Anderson

Andrej

Barker

, et al. (2021) MFEM: a modular finite element methods library. Computers & Mathematics with Applications 81: 42–74. https://doi.org/10.1016/j.camwa.2020.06.009

Arndt

Bangerth

Davydov

, et al. (2021) The deal.II finite element library: design, features, and insights. Computers & Mathematics with Applications 81: 407–422. Development and Application of Open-source Software for Problems with Numerical PDEs. https://doi.org/10.1016/j.camwa.2020.02.022

Bailey

Filippelli

Baker

(2015) Metocean Data Needs Assessment for U.S. Offshore Wind Energy. Technical report. AWS Truepower LLC.

Balay

Gropp

McInnes

, et al. (1997) Efficient management of parallelism in object oriented numerical software libraries. In: Arge

Bruaset

Langtangen

(eds) Modern Software Tools in Scientific Computing. Birkhäuser Press, pp. 163–202.

Balay

Abhyankar

Adams

, et al. (2024) PETSc/TAO Users Manual. Technical Report ANL-21/39 - Revision 3.22. Argonne National Laboratory.

Beji

Battjes

(1994) Numerical simulation of nonlinear wave propagation over a bar. Coastal Engineering 23(1): 1–16. https://doi.org/10.1016/0378-3839(94)90012-4

10.

Bercea

McRae

ATT

Ham

, et al. (2016) A structure-exploiting numbering algorithm for finite elements on extruded meshes, and its performance evaluation in Firedrake. Geoscientific Model Development 9(10): 3803–3815. https://doi.org/10.5194/gmd-9-3803-2016

11.

Bihs

Wang

Pakozdi

, et al. (2020) REEF3D::FNPF A flexible fully nonlinear potential flow solver. Journal of Offshore Mechanics and Arctic Engineering 142(4): 041902. https://doi.org/10.1115/1.4045915

12.

Bingham

Zhang

(2007) On the accuracy of finite-difference solutions for nonlinear water waves. Journal of Engineering Mathematics 58(1): 211–228. https://doi.org/10.1007/s10665-006-9108-4

13.

Bitner-Gregersen

(2015) Joint met-ocean description for design and operations of marine structures. Applied Ocean Research 51: 279–292. https://doi.org/10.1016/j.apor.2015.01.007

14.

Bokhove

Kalogirou

(2016) Variational water wave modelling: from continuum to experiment. In: Bridges

Groves

Nicholls

(eds) Lectures on the Theory of Water Waves, LMS Lecture Note Series. Cambridge University Press, Vol. 426, pp. 226–260. https://doi.org/10.1017/cbo9781316411155.012

15.

Cantwell

Moxey

Comerford

, et al. (2015) Nektar++: an open-source spectral/hp element framework. Computer Physics Communications 192: 205–219. https://doi.org/10.1016/j.cpc.2015.02.008

16.

Chalmers

Karakus

Austin

, et al. (2022) libParanumal: A performance portable high-order finite element library. Release 0.4. 0.

17.

Chang

Tsaur

Huang

(2012) Accurate solution to diffraction around a modified v-shaped breakwater. Coastal Engineering 68: 56–66. https://doi.org/10.1016/j.coastaleng.2012.05.002

18.

Christiansen

Bingham

Engsig-Karup

, et al. (2013) Efficient hybrid-spectral model for fully nonlinear numerical wave tank. In: Proceedings of the 32nd International Conference on Ocean. Offshore and Arctic Engineering.

19.

Dean

Dalrymple

(1991) Water Wave Mechanics for Engineers and Scientists. 2nd edition. World Scientific.

20.

Dommermuth

Yue

DKP

(1987) A high-order spectral method for the study of nonlinear gravity waves. Journal of Fluid Mechanics 184: 267–288. https://doi.org/10.1017/s002211208700288x

21.

Ducrozet

Bingham

Engsig-Karup

, et al. (2010) High-order finite difference solution for 3D nonlinear wave-structure interaction. Journal of Hydrodynamics 22(1): 225–230. https://doi.org/10.1016/s1001-6058(09)60198-0

22.

Engsig-Karup

(2006) Unstructured Nodal DG-FEM Solution of High-Order Boussinesq-Type Equations. PhD thesis, Department of Mechanical Engineering, Technical University of Denmark.

23.

Engsig-Karup

(2014) Analysis of efficient preconditioned defect correction methods for nonlinear water waves. International Journal for Numerical Methods in Fluids 74(10): 749–773. https://doi.org/10.1002/fld.3873

24.

Engsig-Karup

Eskilsson

(2019) Spectral element FNPF simulation of focused wave groups impacting a fixed FPSO-type body. International Journal of Offshore and Polar Engineering 29: 141–148. https://doi.org/10.17736/ijope.2019.jc753

25.

Engsig-Karup

Laskowski

(2021) An efficient p-multigrid spectral element model for fully nonlinear water waves and fixed bodies. International Journal for Numerical Methods in Fluids 93(9): 2823–2841. https://doi.org/10.1002/fld.5011

26.

Engsig-Karup

Hesthaven

Bingham

, et al. (2006) Nodal DG-FEM solutions of high-order boussinesq-type equations. Journal of Engineering Mathematics 56: 351–370. https://doi.org/10.1007/s10665-006-9064-z

27.

Engsig-Karup

Bingham

Lindberg

(2009) An efficient flexible-order model for 3D nonlinear water waves. Journal of Computational Physics 228(6): 2100–2118. https://doi.org/10.1016/j.jcp.2008.11.028

28.

Engsig-Karup

Eskilsson

Bigoni

(2016) A stabilised nodal spectral element method for fully nonlinear water waves. Journal of Computational Physics 318: 1–21. https://doi.org/10.1016/j.jcp.2016.04.060

29.

Engsig-Karup

Monteserin

Eskilsson

(2019) A mixed Eulerian-Lagrangian spectral element method for nonlinear wave interaction with fixed structures. Water Waves 1(2): 315–342. https://doi.org/10.1007/s42286-019-00018-5

30.

Eskilsson

Abedi

Engsig-Karup

(2025) A co-simulation approach to modelling fully nonlinear water waves using the spectral/hp element method. Ocean Engineering 340: 122385. https://doi.org/10.1016/j.oceaneng.2025.122385

31.

Evans

Newman

(1987) Report on the first international workshop on water waves and floating bodies. Journal of Fluid Mechanics 174: 521–528. https://doi.org/10.1017/s0022112087000223

32.

Falgout

Jones

Yang

(2006) The design and implementation of hypre, a library of parallel high performance preconditioners Numerical Solution of Partial Differential Equations on Parallel Computers. Springer, pp. 267–294.

33.

Fenton

(1988) The numerical solution of steady water wave problems. Computers & Geosciences 14(3): 357–368. https://doi.org/10.1016/0098-3004(88)90066-0

34.

Fenton

(1990) Nonlinear wave theories. Sea 9(1): 3–25.

35.

Geuzaine

Remacle

(2009) Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering 79(11): 1309–1331. https://doi.org/10.1002/nme.2579

36.

Ham

Kelly

PHJ

Mitchell

, et al. (2023) Firedrake User Manual. 1st edition. Imperial College London and University of Oxford and Baylor University and University of Washington.

37.

Harris

Dombre

Benoit

, et al. (2022) Nonlinear time-domain wave-structure interaction: a parallel fast integral equation approach. International Journal for Numerical Methods in Fluids 94(2): 188–222. https://doi.org/10.1002/fld.5051

38.

Hestenes

Stiefel

(1952) Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards 49(6): 409–436. https://doi.org/10.6028/jres.049.044

39.

Hesthaven

Kirby

(2008) Filtering in Legendre spectral methods. Mathematics of Computation 77(263): 1425–1452. https://doi.org/10.1090/s0025-5718-08-02110-8

40.

Kärnä

Kramer

Mitchell

, et al. (2018) Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations. Geoscientific Model Development 11(11): 4359–4382. https://doi.org/10.5194/gmd-11-4359-2018

41.

Karniadakis

Sherwin

(2005) Spectral/hp Element Methods for Computational Fluid Dynamics. OUP.

42.

Kirby

Karniadakis

(2003) De-aliasing on non-uniform grids: algorithms and applications. Journal of Computational Physics 191(1): 249–264. https://doi.org/10.1016/s0021-9991(03)00314-0

43.

Kreiss

Oliger

(1972) Comparison of accurate methods for the integration of hyperbolic equations. Tellus 24(3): 199–215. https://doi.org/10.3402/tellusa.v24i3.10634

44.

Kriebel

(1992) Nonlinear wave interaction with a vertical circular cylinder. Part II: wave run-up. Ocean Engineering 19(1): 75–99. https://doi.org/10.1016/0029-8018(92)90048-9

45.

Kurnia

Ducrozet

(2023) NEMOH: open-source boundary element solver for computation of first- and second-order hydrodynamic loads in the frequency domain. Computer Physics Communications 292: 108885. https://doi.org/10.1016/j.cpc.2023.108885

46.

Larsen

Dancy

(1983) Open boundaries in short wave simulations - a new approach. Coastal Engineering 7(3): 285–297. https://doi.org/10.1016/0378-3839(83)90022-4

47.

Le Roux

Lin

Staniforth

(2000) A semi-implicit semi-lagrangian finite-element shallow-water ocean model. Monthly Weather Review 128(5): 1384–1401. https://doi.org/10.1175/1520-0493(2000)128<1384:asislf>2.0.co;2

48.

Lee

Newman

(2006) Wamit User Manual. WAMIT, Inc.

49.

Liu

Dommermuth

Yue

DKP

(1992) A high-order spectral method for nonlinear wave-body interactions. Journal of Fluid Mechanics 245: 115–136. https://doi.org/10.1017/s0022112092000375

50.

Liu

Xue

Yue

DKP

(2001) Computations of fully nonlinear three-dimensional wave-wave and wave-body interactions. Part 2. Nonlinear waves and forces on a body. Journal of Fluid Mechanics 438: 41–66. https://doi.org/10.1017/s0022112001004384

51.

Yan

(2009) QALE-FEM for numerical modelling of non-linear interaction between 3D moored floating bodies and steep waves. International Journal for Numerical Methods in Engineering 78(6): 713–756. https://doi.org/10.1002/nme.2505

52.

MacCamy

Fuchs

(1954) Wave Forces On Piles: A Diffraction Theory. US Beach Erosion Board, Vol. 69.

53.

McRae

ATT

Bercea

Mitchell

, et al. (2016) Automated generation and symbolic manipulation of tensor product finite elements. SIAM Journal on Scientific Computing 38(5): 25–47. https://doi.org/10.1137/15m1021167

54.

Melander

Engsig-Karup

(2024) A p-multigrid hybrid-spectral model for nonlinear water waves. Water Waves 6(3): 547–570. https://doi.org/10.1007/s42286-024-00098-y

55.

Melander

Laskowski

Sherwin

, et al. (2024) A p-multigrid accelerated nodal spectral element method for free-surface incompressible Navier-Stokes model of nonlinear water waves. ArXiv Preprint: 2411.14977.

56.

Minniti

Visbech

Eskilsson

, et al. (2025) A high-order spectral element solver for steady-state free surface flows. ArXiv Preprint: 2512.23648.

57.

Patera

(1984) A spectral element method for fluid dynamics: Laminar flow in a channel expansion. Journal of Computational Physics 54(3): 468–488. https://doi.org/10.1016/0021-9991(84)90128-1

58.

Ransley

Yan

Brown

, et al. (2019) A blind comparative study of focused wave interactions with a fixed FPSO-like structure (CCP-WSI blind test series 1). International Journal of Offshore and Polar Engineering 29(02): 113–127. https://doi.org/10.17736/ijope.2019.jc748

59.

Rathgeber

Ham

Mitchell

, et al. (2016) Firedrake: automating the finite element method by composing abstractions. ACM Transactions on Mathematical Software 43(3): 1–27. https://doi.org/10.1145/2998441

60.

Robertson

Sherwin

(1999) Free-surface flow simulation using hp/spectral elements. Journal of Computational Physics 155(1): 26–53. https://doi.org/10.1006/jcph.1999.6328

61.

Rønquist

Patera

(1987) A Legendre spectral element method for the Stefan problem. International Journal for Numerical Methods in Engineering 24(12): 2273–2299. https://doi.org/10.1002/nme.1620241204

62.

Seixas de Medeiros

Liu

Yue

(2024) A fast high-order boundary element method for nonlinear water waves generation and propagation in large wave basins. Computer Methods in Applied Mechanics and Engineering 432: 117396. https://doi.org/10.1016/j.cma.2024.117396

63.

Shao

Faltinsen

(2014) A harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics. Journal of Computational Physics 274: 312–332. https://doi.org/10.1016/j.jcp.2014.06.021

64.

Thomas

Dwarakish

(2015) Numerical wave modelling – a review. Aquatic Procedia 4: 443–448. https://doi.org/10.1016/j.aqpro.2015.02.059

65.

Tong

Shao

Bingham

, et al. (2024) An adaptive harmonic polynomial cell method for three-dimensional fully nonlinear wave-structure interaction with immersed boundaries. Physics of Fluids 36(3): 032118. https://doi.org/10.1063/5.0190037

66.

Visbech

Bingham

Eskilsson

, et al. (2024a) A high-order accurate spectral element based time-domain simulation of a model-scale floating offshore wind turbine. International Journal of Offshore and Polar Engineering 34: 254–262. https://doi.org/10.17736/ijope.2024.sv17

67.

Visbech

Engsig-Karup

Bingham

(2024b) Solving the complete pseudo-impulsive radiation and diffraction problem using a spectral element method. Computer Methods in Applied Mechanics and Engineering 423: 116871. https://doi.org/10.1016/j.cma.2024.116871

68.

Visbech

Engsig-Karup

Bingham

, et al. (2026c) Recent progress on modeling nonlinear wave propagation and wave-structure interaction using a high-order shifted boundary method: capabilities, challenges, and perspectives. In: Proceedings of the 41th International Workshop on Water Waves and Floating Bodies.

69.

Willberg

Duczek

Perez

, et al. (2012) Comparison of different higher order finite element schemes for the simulation of Lamb waves. Computer Methods in Applied Mechanics and Engineering 241: 246–261.

70.

Eatock Taylor

(1994) Finite element analysis of two-dimensional non-linear transient water waves. Applied Ocean Research 16(6): 363–372. https://doi.org/10.1016/0141-1187(94)00029-8

71.

(2004) Simulation of nonlinear interactions between waves and floating bodies through a finite-element-based numerical tank. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences 460(2050): 2797–2817. https://doi.org/10.1098/rspa.2004.1302

72.

Cantwell

Monteserin

, et al. (2018) Spectral/hp element methods: recent developments, applications, and perspectives. Journal of Hydrodynamics 30: 1–22. https://doi.org/10.1007/s42241-018-0001-1

73.

Xue

Liu

, et al. (2001) Computations of fully nonlinear three-dimensional wave-wave and wave-body interactions. Part 1. Dynamics of steep three-dimensional waves. Journal of Fluid Mechanics 438: 11–39. https://doi.org/10.1017/s0022112001004396

74.

Yan

Liu

(2011) An efficient high-order boundary element method for nonlinear wave-wave and wave-body interactions. Journal of Computational Physics 230(2): 402–424. https://doi.org/10.1016/j.jcp.2010.09.029

75.

Zakharov

(1968) Stability of periodic waves of finite amplitude on the surface of a deep fluid. Journal of Applied Mechanics and Technical Physics 9: 190–194. https://doi.org/10.1007/bf00913182