category: literaturenote citekey: visbechfnpfsemparallelspectral2025 title: "FNPF-SEM: A parallel spectral element model in Firedrake for fully nonlinear water wave simulations" authors: "Visbech, Jens; Melander, Anders; Engsig-Karup, Allan Peter" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.09435v1" zotero_key: ERQHCZDE zotero_storage: 8W9GXWIN collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Visbech, Jens"

status: converted

FNPF-SEM: A PARALLEL SPECTRAL ELEMENT MODEL IN FIREDRAKE FOR FULLY NONLINEAR WATER WAVE SIMULATIONS

AN ARXIV PREPRINT

Jens Visbech

Department of Applied Mathematics and Computer Science Technical University of Denmark Kongens Lyngby, 2800, Denmark jvis@dtu.dk

Anders Melander

Department of Applied Mathematics and Computer Science Technical University of Denmark Kongens Lyngby, 2800, Denmark adame@dtu.dk

Allan Peter Engsig-Karup

Department of Applied Mathematics and Computer Science Technical University of Denmark Kongens Lyngby, 2800, Denmark apek@dtu.dk

June 12, 2025

ABSTRACT

We present a new parallel spectral element solver, FNPF-SEM, for simulating linear and fully nonlinear potential flow-based water waves and their interaction with offshore structures. The tool is designed as a general-purpose wave model for offshore engineering applications. Built within the open-source framework Firedrake, the new FNPF-SEM model is designed as a computational tool capable of capturing both linear and nonlinear wave phenomena with high accuracy and efficiency, with support for high-order (spectral) finite elements. Additionally, Firedrake provides native support for MPI-based parallelism, allowing for efficient multi-CPU distributed computations needed for largescale simulations. We demonstrate the capabilities of the high-order spectral element model through h- and p-convergence studies, and weak and strong scaling tests. Validation is performed against analytical solutions and experimental data for several benchmark cases, including nonlinear highorder harmonic generation and linear and nonlinear 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 Spectral element method · linear and nonlinear free surface potential flow · 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 considering complex geometry due to the flow domain and included marine or offshore structures. 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 [6], wave energy devices, and breakwaters. Other cases include metocean studies such as assessments of ocean conditions, [12]. Many offshore engineering applications have to consider modeling of nonlinear water waves in large regions and within large time scales, requiring effective parallel implementations of such wave models to accurately and efficiently compute real-life scenarios within acceptable computational times.

Physical accuracy is typically associated with the chosen mathematical model and how it represents 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 flow applications can be approximated by fully nonlinear potential flow (FNPF), assuming the fluid is incompressible, inviscid, and irrotational. 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 complexity. One desirable feature of the FNPF formulation is the ability to simulate both offshore and nearshore environments, since wave interactions with the bathymetry are accurately represented. 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 [50]. 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 the FNPF and LPF formulations have been a major research topic for decades. 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 [32, 40, 42, 44, 53, 61, 62]. Another option are high-order spectral methods, which are very efficient when modeling nonlinear wave propagation and submerged wave-structure interactions, [17, 18, 43, 47]. The FNPF and LPF equations have also been solved using finite difference schemes, [2, 3, 11, 19, 21], and using harmonic polynomials, [54, 55]. Using the classical second-order accurate finite element method, we have the work by [59] and many more.

Spectral element methods (SEMs), initially introduced by [49], 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 and highly accurate numerical solutions due to the use of 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 traditional numerical approaches' limitations, e.g., difficulty in modeling complex geometry and/or high dispersive and diffusive errors. For further details and applications of the SEM, see the review by [60] and the textbook by [36].

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++, [14], the arbitrary high-order finite element framework MFEM, [4], the deal.II finite element framework, [5], or libParanumal, [15]. Another computational framework that allows for high-order basis functions is Firedrake [31]. Firedrake is an open-source framework that leverages itself around automating the required steps to solve partial differential equations using finite elements. One of Firedrake's main features is the unified form language, [1], which allows for the symbolic expression of, e.g., any weak variational formulations. Moreover, Firedrake uses state-of-the-art PETSc solvers [7, 8], and supports structured and unstructured meshes using the most common types of elements and finite element spaces. Lastly, Firedrake allows for a fully parallelized computational setup through MPI and can run on large many-core hardware 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 [23], who proposed a solution to the intrinsic mesh-induced instability associated with finite element discretizations of free surface flows, as previously identified in [52]. The SEM implementations at DTU are based primarily on the FNPF formulations, [22, 23, 24, 25, 58], and the LPF approximations, [56, 57]. This combined FNPF-SEM approach offers a numerical approach that can handle complex geometries and simulate over long time periods accurately by using the high-order polynomial basis. The methodology has been extended to encompass the full incompressible Navier–Stokes equations in [48] with a free surface. Achieving low dispersive and diffusive numerical errors that stem from the numerical discretization procedure is a key aspect for long-time wave simulations. High-order discretization methods, such as SEM, are very efficient compared to low-order methods [23, 38]. In addition, [26] 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 and scalable features provided by modern open-source computational frameworks. Moreover, efficient water wave solutions to the Benney-Luke-type equations stated using variational principles are also available via Firedrake, [13]. Furthermore, the Thetis project, an unstructured grid coastal ocean model, was built using the Firedrake framework, [35].

1.1 Paper contribution

The main aim is to verify and validate 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, allowing for an efficient parallel implementation of the FNPF equations to support large- and full-scale wave problems over long time horizons, ultimately to be used within 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 validated in various test cases for FNPF and LPF wave propagation and wave structures 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 approaches, highlighting the chosen spatial and temporal discretization approaches, including key aspects such as wave initialization, generation, and absorption. The results are shown in Section 4, which presents both numerical verification studies (including convergence and scaling analyses) and validation cases involving 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 is also included for supplementary material in Section 6.

2 Mathematical Model

Considering a fluid domain, Ω ∈ 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: Γ = S i Γ i for i ∈ {η, b, w, p}. Here, Γ η is the top free surface boundary defined at z = η in the FNPF setting, where η = η(x, y, t) : Γ^η × T 7→ R is the free surface elevation defined from z = 0 in the time-domain T : t ≥ 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 7→ 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) : Ω × T 7→ R, can be used to express the fluid velocities, u = (u, v, w) = (u(x, t), v(x, t), w(x, t)) : Ω × T 7→ 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

$$\nabla^{2} \phi = 0, \quad \text{in} \quad \Omega,$$

$$\phi = \phi_{\eta}, \quad \text{on} \quad \Gamma^{\eta},$$

$$\nabla \phi \cdot \mathbf{n} = 0, \quad \text{on} \quad {\Gamma^{b}, \Gamma^{w}},$$

$$\phi|{\Gamma^{p}} = \phi|{\Gamma^{p}}, \quad \text{on} \quad \Gamma^{p},$$

$$(1)$$

where n = (nx, ny, nz) = (nx(x), ny(x), nz(x)) : Γ 7→ 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) : Γ^η × T 7→ R. This notation becomes useful when we consider the FNPF time governing kinematic and dynamic free surface conditions on Zakharov form, [63]. On Γ ^η × T , we have

$$\partialt \eta = -\widetilde{\nabla} \eta \cdot \widetilde{\nabla} \phi{\eta} + w_{\eta} (1 + \widetilde{\nabla} \eta \cdot \widetilde{\nabla} \eta), \tag{2}$$

$$\partialt \phi{\eta} = -g\eta - \frac{1}{2} \widetilde{\mathbf{\nabla}} \phi{\eta} \cdot \widetilde{\mathbf{\nabla}} \phi{\eta} + \frac{1}{2} w_{\eta}^2 (1 + \widetilde{\mathbf{\nabla}} \eta \cdot \widetilde{\mathbf{\nabla}} \eta). \tag{3}$$

Here ^∇^e = (∂x, ∂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. Remark: The classical linearized version of (1) through (3) is considered as well; however, these governing equations are omitted for brevity.

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, [31]. 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, [1], allowing users to express weak variational formulations in a compact and 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, [7], 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 mesh tessellations.

3.1 Spatial discretization

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

$$\mathcal{D} \simeq \mathcal{T}{\mathcal{D}} = \bigcup{n=1}^{N_{\mathcal{D}}} \mathcal{D}_n, \quad \text{for} \quad \mathcal{D} \in {\Omega, \Gamma^{\eta}},$$

where N^D is an arbitrary number of elements, such that Dⁿ is the n'th element in the tesselation TD. The element types and meshes are highlighted in the following section. On the two tesselations, we introduce C (D) function spaces of piecewise polynomial functions as:

$$\mathcal{V}{\mathcal{D}} = { \nu{\mathcal{D}} \in C^{0}(\mathcal{D}); \forall n \in {1, ..., N{\mathcal{D}}}, \nu{\mathcal{D}}|{\mathcal{D}{n}} \in \mathbb{P}^{p} },$$

where P ^p denotes the polynomial space of nodal Lagrange polynomials of degree at most p.

3.1.1 Weak variational formulations

Next, we state the weak forms of (1) and (2) through (3), respectively. The linear variational problem for the Laplace problem is to find ϕ ∈ V^Ω such that ∀ν^Ω ∈ VΩ, with

$$a[\phi; \nu{\Omega}] = b[\nu{\Omega}], \tag{4}$$

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

$$a[\phi; \nu{\Omega}] = \langle \nabla \phi, \nabla \nu{\Omega} \rangle_{\Omega},$$

$b[\nu_{\Omega}] = 0,$

with ⟨a, b⟩^S = R 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 {∂tη, ∂tϕη} ∈ VΓ^η such that ∀νΓ^η ∈ VΓ^η , we solve

$$F_k[\partialt \eta; \nu{\Gamma^{\eta}}] = 0$$ , and $F_d[\partialt \phi{\eta}; \nu_{\Gamma^{\eta}}] = 0$ , (5)

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

$$F_k[\partialt \eta; \nu{\Gamma^{\eta}}] = \left\langle \partialt \eta, \nu{\Gamma^{\eta}} \right\rangle{\Gamma^{\eta}} + \left\langle \widetilde{\nabla} \eta \cdot \widetilde{\nabla} \phi{\eta} - w{\eta} (1 + \widetilde{\nabla} \eta \cdot \widetilde{\nabla} \eta), \nu{\Gamma^{\eta}} \right\rangle_{\Gamma^{\eta}},$$

and

$$F_d[\partialt \phi{\eta}; \nu_{\Gamma^{\eta}}] = \left\langle \partialt \phi{\eta}, \nu{\Gamma^{\eta}} \right\rangle{\Gamma^{\eta}} + \left\langle g \eta - \frac{1}{2} \widetilde{\boldsymbol{\nabla}} \phi{\eta} \cdot \widetilde{\boldsymbol{\nabla}} \phi{\eta} - \frac{1}{2} w{\eta}^2 (1 + \widetilde{\boldsymbol{\nabla}} \eta \cdot \widetilde{\boldsymbol{\nabla}} \eta), \nu{\Gamma^{\eta}} \right\rangle_{\Gamma^{\eta}}.$$

3.1.2 Initial meshes, extrusion, and mesh updating

Firedrake has native support for extruded meshes, [10, 46]. 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, [30]. 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. Third, in the FNPF setting, a vertical translation is applied to the mesh along each column of points to track the evolving free surface curvature. This ensures an accurate geometric representation of the time-dependent fluid domain, Ω, as 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. Remark: A similar approach is applied 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 being done to include this.

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^Ω 7→ 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.

Figure 2: Point distribution (indicated by different colors) for a finite FNPF wave tank distributed over N^c = 5 cores on a mesh with (Nx, Ny, Nz) = (20, 4, 4) elements of order p = 5.

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 = CCFL∆xminu −1 max. Here, CCFL = 0.95, ∆xmin denotes the minimum spacing between two grid points on Γ η , and umax represents the maximum velocity encountered in the simulation.

Some form of stabilization is typically required to model strongly nonlinear wave phenomena accurately. In [23], the instability issue initially highlighted by [52] 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, [37], and a mild modal (spectral) filter applied to the free surface quantities, [34]. 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, [23]. 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 for the parallel application of the filter across all free surface elements on the interior point values. Specifically, the kernel operates element-wise, applying the local modal filter as defined in [23]: Suppose gn(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ⁿ = (xn, yn). We map the solution to the modal space as

$$\widehat{g}_n(f) = \mathcal{V}^{-1} g_n(\boldsymbol{x_n}),$$

where V ∈ R ⁿp×n^p is a conforming Vandermonde matrix defined as ^Vij ⁼ ^ψb^j (ri), where ^rⁱ = (rⁱ , si), is coordinates in a reference domain, ^R ⁼ {^r = (r, s)|(r, s) ≥ −1; ^r ⁺ ^s ^≤ ⁰} and ^ψb(r) is a modal basis function. See [23] for more details. Using a filter matrix, F ∈ R ⁿp×n^p , filtering can be applied to the highest modes

$$\widehat{g}_{n,filtered}(f) = \mathcal{F} \widehat{g}_n(f).$$

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

$$g_{n,filtered}(\boldsymbol{xn}) = \mathcal{V} \ \widehat{g}{n,filtered}(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, [28]. 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, [29]. 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, [41], 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 + Cr f{\text{true}},$$

where f ∈ {η, ϕη} and C^r = Cr(x, y) : Γ^η 7→ 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, ftrue 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 a duration of 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 = tstart.

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

Step 4.1 – Mesh-update: Update T^Ω according to η.

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 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 related to the verification and validation of the model. All simulations are performed on the HPC system at the DTU Computing Center at the Technical University of Denmark using the hardware configurations listed in Table 1.

| Total nodes, Nn, available | 15 | | |----------------------------|-----------------------|--| | Total cores, Nc, available | 300 | | | Cores per node, Nc/Nn | 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 [51] by using a CG solver with the HYPRE BoomerAMG algebraic multigrid preconditioner, [27, 33], 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 solver with a ILU as preconditioner and a relative and absolute tolerance of 10⁻⁵ and 10⁻¹⁵, respectively. For the purpose of 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 [20].

4.1 Model verification

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

4.1.1 Convergence hp 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η|∞, 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 hmax 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, [23], leading to an expected convergence rate of O(h p max). This expected behavior is clearly observed in the results. Furthermore, increasing the nonlinearity of the wave is shown to increase the complexity of the problem, 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 A, in Figure 15, where spectral, i.e., exponential, convergence is observed. The test is conducted on six different meshes with increasing resolution. Specifically, Mesh 1 has hmax = 0.25, Mesh 2 has hmax = 0.125, Mesh 3 has hmax = 0.0625, etc. As observed previously, increasing nonlinearity results in a more challenging problem, requiring higher resolution to achieve comparable accuracy. Remark: Similar analyses have been performed for the LPF model, which exhibits improved convergence behavior, which is an expected outcome for a linearized model. These 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 per time stage for 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, from comparing the results of Mesh 1 and 2, it is clear that increasing the problem size allows for scalable behavior on a higher amount of cores due to LaplaceSolve being the dominant cost factor for longer.

Figure 4: Timings for the three dominant 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 due to the constant cost of the remaining routines. 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 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 (4) and the nonlinear variational problem in (5), which reduces to a simple function evaluation in the linearized setting.

Figure 5: Percentage distribution of computational cost across the three dominant routines for Mesh 1 (blue) and Mesh 2 (red), using p = 4.

Based on this performance profiling of the code, it is confirmed that the Laplace solver constitutes 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, is has been for conciseness.

4.1.3 Parallel performance: Strong scaling

To evaluate how run time performance is affected by increasing the number of cores for a fixed problem size (i.e., 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, Nn, where each node comprises 20 cores. As a metric for strong scaling efficiency, we define γ^s as

$$\gamma_s = \frac{T^b Nn^b}{T{N_n} N_n},\tag{6}$$

where T b and N^b ⁿ denote the runtime and node count of the baseline simulation, respectively, and TNⁿ is the runtime using Nⁿ nodes. In this study, the baseline corresponds to a single node configuration (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ⁿ = 15 nodes, equivalent to N^c = 300 cores, using a nonlinear periodic domain with kh = 2π.

Table 3: Strong scaling meshes with Nz = 4.
Mesh	NΩ	NΓη	p	DoF(TΩ)	DoF(TΓη )
			3	479,232	36,864
1	32,768	8,192	4	1,114,112	65,536
			5	2,150,400	102,400
	131,072	32,768	3	1,916,928	147,456
2			4	4,456,448	262,144
			5	8,601,600	409,600
			3	3,833,037	294,849
	3 262,088	65,522	4	8,910,992	524,176
			5	17,199,525	819,025

Figure 6 shows the strong scaling efficiency using (6) for the three meshes. As expected, both increased mesh refinement and higher polynomial order improve the parallel efficiency of the model. 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ⁿ = 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 performed 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ⁿ = {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 = NcL, 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

$$\gammaw = \frac{T^b}{T{N_n}}. (7)$$

Mesh	NΩ	Nx	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

Table 4: Weak scaling meshes with N^y = N^z = 4 and polynomial order p = 5. All reported quantities are per core, Nc.

Figure 8 presents the run times for the various mesh configurations. As expected, increasing resolution leads to longer run times. 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 (7). 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 per stage on six meshes with polynomial order p = 5.

Figure 9: Weak scaling efficiency on six meshes with polynomial order p = 5.

4.2 Model validation

For all the validation 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 against the physical experimental laboratory data by Beji and Battjes, [9], 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 make it fit 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). Absorption zone (30 ≤ x ≤ 38). The width of the domain is 0 ≤ y ≤ 1. All units are [m]. The depth before and after the bar is h = 0.4 [m] and h = 0.1 [m] on top of 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 7,600 elements, which is extruded with N^z = 4 elements in the vertical to give 22,800 prismatic elements for the representation of Ω. Ultimately, this gives a problem size 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.

Figure 10: Contour plot of the free surface elevation at t = 25T for the bar test.

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

A contour plot of the free surface elevation, η, at the final time is shown in Figure 10. Here, the high-order harmonic can be seen to take place 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πfnt, where fⁿ = 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.

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 surfacepiercing 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 B, Figure 16, the unstructured mesh used to obtain the linear solution can be seen.

	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

Table 5: Domain, mesh, and wave parameters for vertical cylinder simulations.

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

$$\eta_m(x_i, y_i) = \max_t \eta(x_i, y_i, t), \quad (x_i, y_i) \in \Gamma^{\text{cylinder}}$$

For the linear test case (TC 1), the numerical results, denoted LPF-1, are compared with the analytical predictions of linear diffraction theory by [45], 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 [39], denoted K92-2, which were digitized for comparison. We also include equivalent linear results for this configuration, denoted MF54-2.

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

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 [19]. 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 validation test case, we consider wave interactions with a breakwater structure, following the setup described in [16], 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. 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. While the original study by [16] treated the breakwaters as infinitely thin, we adopt a more physically realistic representation by assigning them a finite thickness of l = 1 [m].

Figure 13: Free surface mesh for θ = 0◦ around the breakwater.

The computational domain spans the region −L² − 2b ≤ (x, y) ≤ 2b + L2, completely enclosing the target region defined as −2b ≤ x ≤ 2b and −L² − 2b ≤ y ≤ 2b + L2. Here, L² = 260 [m] ≈ 2L, 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 + L2. 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 ten 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). The free surface mesh can be seen in Figure 13.

The wave input parameters are taken from [16], with a non-dimensional wavenumber kh = 0.65 [-], corresponding to L = 129.5 [m]. The linear – and arbitrary – wave height is set to H = 1.295 [m], and the simulations are run for 10 wave periods. For each incident angle, θ = {0 ◦ , 30◦ , 60◦}, we compute the time-independent free surface elevation variation, defined at each grid point (xⁱ , yi) as

$$\eta_v = \eta_v(x_i, y_i) = \max_t \eta(x_i, y_i, t) - \min_t \eta(x_i, y_i, t).$$

The corresponding non-dimensional quantity ηv/H is visualized as contour plots in Appendix C (Figures 17, 18, and 19). 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 [16] (C12-θ). Despite small differences arising from the use of finite-thickness structures, the discretizations of the domain, and minor wave reflections at the boundaries, the numerical solutions show good qualitative agreement with the reference results. This highlights that the incident angle impacts the wave elevation in the vicinity of the breakwater structure, e.g., which may be of interest in 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 [16] (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 largescale computations efficiently. We detailed the mathematical formulation, including the governing equations, domain definitions, and numerical discretization. Parallel performance was assessed through code performance profiling for strong- and weak-scaling studies, which demonstrated 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 validated 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 interaction 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 consider areas of irregular coastlines, exploiting the advantage of unstructured meshes.

Acknowledgments

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.

Declarations

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

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".

References

[1] Alnæs MS, Logg A, Ølgaard KB, Rognes ME and Wells GN (2014) Unified form language: A domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software 40(2).
[2] Amini-Afshar M and Bingham HB (2017) Solving the linearized forward-speed radiation problem using a high-order finite difference method on overlapping grids. Applied Ocean Research 69: 220–244.
[3] Amini-Afshar M and Bingham HB (2018) Pseudo-impulsive solutions of the forward-speed diffraction problem using a high-order finite-difference method. Applied Ocean Research 80: 197–219.
[4] Anderson R, Andrej J, Barker A, Bramwell J, Camier JS, Cerveny J, Dobrev V, Dudouit Y, Fisher A, Kolev T, Pazner W, Stowell M, Tomov V, Akkerman I, Dahm J, Medina D and Zampini S (2021) MFEM: A modular finite element methods library. Computers & Mathematics with Applications 81: 42–74.
[5] Arndt D, Bangerth W, Davydov D, Heister T, Heltai L, Kronbichler M, Maier M, Pelteret JP, Turcksin B and Wells D (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.
[6] Bailey BH, Filippelli M and Baker M (2015) Metocean data needs assessment for u.s. offshore wind energy. Technical report, AWS Truepower LLC., Albany, NY (United States).
[7] Balay S, Abhyankar S, Adams MF, Benson S, Brown J, Brune P, Buschelman K, Constantinescu E, Dalcin L, Dener A, Eijkhout V, Faibussowitsch J, Gropp WD, Hapla V, Isaac T, Jolivet P, Karpeev D, Kaushik D, Knepley MG, Kong F, Kruger S, May DA, McInnes LC, Mills RT, Mitchell L, Munson T, Roman JE, Rupp K, Sanan P, Sarich J, Smith BF, Suh H, Zampini S, Zhang H, Zhang H and Zhang J (2024) PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.22, Argonne National Laboratory.
[8] Balay S, Gropp WD, McInnes LC and Smith BF (1997) Efficient management of parallelism in object oriented numerical software libraries. In: Arge E, Bruaset AM and Langtangen HP (eds.) Modern Software Tools in Scientific Computing. Birkhäuser Press, pp. 163–202.
[9] Beji S and Battjes J (1994) Numerical simulation of nonlinear wave propagation over a bar. Coastal Engineering 23(1): 1–16.
[10] Bercea GT, McRae ATT, Ham DA, Mitchell L, Rathgeber F, Nardi L, Luporini F and Kelly PHJ (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.
[11] Bingham HB and Zhang H (2007) On the accuracy of finite-difference solutions for nonlinear water waves. Journal of Engineering Mathematics 58(1): 211–228.
[12] Bitner-Gregersen EM (2015) Joint met-ocean description for design and operations of marine structures. Applied Ocean Research 51: 279–292.
[13] Bokhove O and Kalogirou A (2016) Variational water wave modelling: from continuum to experiment. In: Bridges TJ, Groves MD and Nicholls DP (eds.) Lectures on the theory of water waves, LMS Lecture Note Series, volume 426. Cambridge University Press, pp. 226–260.
[14] Cantwell C, Moxey D, Comerford A, Bolis A, Rocco G, Mengaldo G, De Grazia D, Yakovlev S, Lombard JE, Ekelschot D, Jordi B, Xu H, Mohamied Y, Eskilsson C, Nelson B, Vos P, Biotto C, Kirby R and Sherwin S (2015) Nektar++: An open-source spectral/hp element framework. Computer Physics Communications 192: 205–219.
[15] Chalmers N, Karakus A, Austin AP, Swirydowicz K and Warburton T (2022) libParanumal: a performance portable high-order finite element library. Release 0.5.0.
[16] Chang KH, Tsaur DH and Huang LH (2012) Accurate solution to diffraction around a modified v-shaped breakwater. Coastal Engineering 68: 56–66.
[17] Christiansen TB, Bingham HB, Engsig-Karup AP, Ducrozet G and Ferrant P (2013) Efficient hybrid-spectral model for fully nonlinear numerical wave tank. In: Proceedings of the 32nd International Conference on Ocean, Offshore and Arctic Engineering.
[18] Dommermuth DG and Yue DKP (1987) A high-order spectral method for the study of nonlinear gravity waves. Journal of Fluid Mechanics 184: 267–288.
[19] Ducrozet G, Bingham HB, Engsig-Karup AP and Ferrant P (2010) High-order finite difference solution for 3D nonlinear wave-structure interaction. Journal of Hydrodynamics 22(1): 225–230.
[20] Engsig-Karup AP (2014) Analysis of efficient preconditioned defect correction methods for nonlinear water waves. International Journal for Numerical Methods in Fluids 74(10): 749–773.
[21] Engsig-Karup AP, Bingham HB and Lindberg O (2009) An efficient flexible-order model for 3D nonlinear water waves. Journal of Computational Physics 228(6): 2100–2118.
[22] Engsig-Karup AP and Eskilsson C (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.
[23] Engsig-Karup AP, Eskilsson C and Bigoni D (2016) A stabilised nodal spectral element method for fully nonlinear water waves. Journal of Computational Physics 318: 1–21.
[24] Engsig-Karup AP and Laskowski WL (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.
[25] Engsig-Karup AP, Monteserin C and Eskilsson C (2019) A Mixed Eulerian–Lagrangian Spectral Element Method for Nonlinear Wave Interaction with Fixed Structures. Water Waves 1(2): 315–342.
[26] Eskilsson C and Engsig-Karup AP (2024) Water wave simulations using fully nonlinear potential flow: Spectral/hp element models implemented in nektar++. In: Proceedings of The 26th Numerical Towing Tank Symposium.
[27] Falgout RD, Jones JE and Yang UM (2006) The design and implementation of hypre, a library of parallel high performance preconditioners. In: Numerical solution of partial differential equations on parallel computers. Springer, pp. 267–294.
[28] Fenton J (1988) The numerical solution of steady water wave problems. Computers & Geosciences 14(3): 357–368.
[29] Fenton JD (1990) Nonlinear wave theories. the Sea 9(1): 3–25.
[30] Geuzaine C and Remacle JF (2009) Gmsh: A 3-D finite element mesh generator with built-in pre- and postprocessing facilities. International Journal for Numerical Methods in Engineering 79(11): 1309–1331.
[31] Ham DA, Kelly PHJ, Mitchell L, Cotter CJ, Kirby RC, Sagiyama K, Bouziani N, Vorderwuelbecke S, Gregory TJ, Betteridge J, Shapero DR, Nixon-Hill RW, Ward CJ, Farrell PE, Brubeck PD, Marsden I, Gibson TH, Homolya M, Sun T, McRae ATT, Luporini F, Gregory A, Lange M, Funke SW, Rathgeber F, Bercea GT and Markall GR (2023) Firedrake User Manual. Imperial College London and University of Oxford and Baylor University and University of Washington, first edition edition.
[32] Harris JC, Dombre E, Benoit M, Grilli ST and Kuznetsov KI (2022) Nonlinear time-domain wave-structure interaction: A parallel fast integral equation approach. International Journal for Numerical Methods in Fluids 94(2): 188–222.
[33] Hestenes MR, Stiefel E et al. (1952) Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards 49(6): 409–436.
[34] Hesthaven JS and Kirby RM (2008) Filtering in legendre spectral methods. Mathematics of Computation 77(263): 1425–1452.
[35] Kärnä T, Kramer SC, Mitchell L, Ham DA, Piggott MD and Baptista AM (2018) Thetis coastal ocean model: discontinuous galerkin discretization for the three-dimensional hydrostatic equations. Geoscientific Model Development 11(11): 4359–4382.
[36] Karniadakis G and Sherwin S (2005) Spectral/hp element methods for computational fluid dynamics. OUP Oxford.
[37] Kirby RM and Karniadakis GE (2003) De-aliasing on non-uniform grids: algorithms and applications. Journal of Computational Physics 191(1): 249–264.
[38] Kreiss HO and Oliger J (1972) Comparison of accurate methods for the integration of hyperbolic equations. Tellus 24(3): 199–215.
[39] Kriebel D (1992) Nonlinear wave interaction with a vertical circular cylinder. part ii: Wave run-up. Ocean Engineering 19(1): 75–99.
[40] Kurnia R and Ducrozet G (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.
[41] Larsen J and Dancy H (1983) Open boundaries in short wave simulations — a new approach. Coastal Engineering 7(3): 285–297.
[42] Lee CH and Newman JN (2006) Wamit user manual. WAMIT, Inc 42.
[43] Liu Y, Dommermuth DG and Yue DKP (1992) A high-order spectral method for nonlinear wave–body interactions. Journal of Fluid Mechanics 245: 115–136.
[44] Liu Y, Xue M and 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.
[45] MacCamy RC and Fuchs RA (1954) Wave forces on piles: a diffraction theory. 69. US Beach Erosion Board.
[46] McRae ATT, Bercea GT, Mitchell L, Ham DA and Cotter CJ (2016) Automated generation and symbolic manipulation of tensor product finite elements. SIAM Journal on Scientific Computing 38(5): S25–S47.
[47] Melander A and Engsig-Karup AP (2024) A p-multigrid hybrid-spectral model for nonlinear water waves. Water Waves 6(3): 547–570.
[48] Melander A, Laskowski W, Sherwin SJ and Engsig-Karup AP (2024) A p-multigrid accelerated nodal spectral element method for free-surface incompressible navier-stokes model of nonlinear water waves. arXiv Preprint: 2411.14977.
[49] Patera AT (1984) A spectral element method for fluid dynamics: Laminar flow in a channel expansion. Journal of Computational Physics 54(3): 468–488.
[50] Ransley E, Yan S, Brown SA, Mai T, Graham D, Ma Q, Musiedlak PH, Engsig-Karup AP, Eskilsson C, Li Q 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.
[51] Rathgeber F, Ham DA, Mitchell L, Lange M, Luporini F, Mcrae ATT, Bercea GT, Markall GR and Kelly PHJ (2016) Firedrake: Automating the finite element method by composing abstractions. ACM Transactions on Mathematical Software 43(3).
[52] Robertson I and Sherwin S (1999) Free-Surface Flow Simulation Using hp/Spectral Elements. Journal of Computational Physics 155(1): 26–53.
[53] Seixas de Medeiros J, Liu Y and Yue DK (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.
[54] Shao YL and Faltinsen OM (2014) A harmonic polynomial cell (hpc) method for 3d laplace equation with application in marine hydrodynamics. Journal of Computational Physics 274: 312–332.
[55] Tong C, Shao Y, Bingham HB and Hanssen FCW (2024) An adaptive harmonic polynomial cell method for three-dimensional fully nonlinear wave-structure interaction with immersed boundaries. Physics of Fluids 36(3): 032118.
[56] Visbech J, Bingham HB, Eskilsson C, Palm J and Engsig-Karup AP (2024) 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.
[57] Visbech J, Engsig-Karup AP and Bingham HB (2024) Solving the complete pseudo-impulsive radiation and diffraction problem using a spectral element method. Computer Methods in Applied Mechanics and Engineering 423: 116871.
[58] Visbech J, Engsig-Karup AP, Bingham HB, Amini-Afshar M and Ricchiuto M (2024) A high-order shifted boundary method for water waves and floating bodies. In: Proceedings of The 39th International Workshop on Water Waves and Floating Bodies.
[59] Wu G and Eatock Taylor R (1994) Finite element analysis of two-dimensional non-linear transient water waves. Applied Ocean Research 16(6): 363–372.
[60] Xu H, Cantwell C, Monteserin C, Eskilsson C, Engsig-Karup AP and Sherwin S (2018) Spectral/hp element methods: Recent developments, applications, and perspectives. Journal of Hydrodynamics 30: 1–22.
[61] Xue M, Xu H, Liu Y and Yue DKP (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.
[62] Yan H and Liu Y (2011) An efficient high-order boundary element method for nonlinear wave–wave and wave-body interactions. Journal of Computational Physics 230(2): 402–424.
[63] Zakharov VE (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.

6 Appendix

6.1 Appendix A: p-convergence study

Figure 15: Convergence study for p for nonlinear waves in all combinations between kh = {1, 3, 6} and ε/εmax = {10%, 50%, 90%}.

6.2 Appendix B: Vertical cylinder mesh

Figure 16: Mesh used for the linear simulation as described in Table 5. Zoomed in around the cylinder.

6.3 Appendix C: Breakwater solutions

Figure 17: Breakwater solution (normalized free surface variation) at θ = 0◦ angle of attack.

Figure 18: Breakwater solution (normalized free surface variation) at θ = 30◦ angle of attack.

Figure 19: Breakwater solution (normalized free surface variation) at θ = 60◦ angle of attack.

visbechfnpfsemparallelspectral2025.md 71 KB Historia Raaka