A Massively Parallel Hybrid Finite Volume/Finite Element Scheme for Computational Fluid Dynamics

Abstract

In this paper, we propose a novel family of semi-implicit hybrid finite volume/finite element schemes for computational fluid dynamics (CFD), in particular for the approximate solution of the incompressible and compressible Navier-Stokes equations, as well as for the shallow water equations on staggered unstructured meshes in two and three space dimensions. The key features of the method are the use of an edge-based/face-based staggered dual mesh for the discretization of the nonlinear convective terms at the aid of explicit high resolution Godunov-type finite volume schemes, while pressure terms are discretized implicitly using classical continuous Lagrange finite elements on the primal simplex mesh. The resulting pressure system is symmetric positive definite and can thus be very efficiently solved at the aid of classical Krylov subspace methods, such as a matrix-free conjugate gradient method. For the compressible Navier-Stokes equations, the schemes are by construction asymptotic preserving in the low Mach number limit of the equations, hence a consistent hybrid FV/FE method for the incompressible equations is retrieved. All parts of the algorithm can be efficiently parallelized, i.e., the explicit finite volume step as well as the matrix-vector product in the implicit pressure solver. Concerning parallel implementation, we employ the Message-Passing Interface (MPI) standard in combination with spatial domain decomposition based on the free software package METIS. To show the versatility of the proposed schemes, we present a wide range of applications, starting from environmental and geophysical flows, such as dambreak problems and natural convection, over direct numerical simulations of turbulent incompressible flows to high Mach number compressible flows with shock waves. An excellent agreement with exact analytical, numerical or experimental reference solutions is achieved in all cases. Most of the simulations are run with millions of degrees of freedom on thousands of CPU cores. We show strong scaling results for the hybrid FV/FE scheme applied to the 3D incompressible Navier-Stokes equations, using millions of degrees of freedom and up to 4096 CPU cores. The largest simulation shown in this paper is the well-known 3D Taylor-Green vortex benchmark run on 671 million tetrahedral elements on 32,768 CPU cores, showing clearly the suitability of the presented algorithm for the solution of large CFD problems on modern massively parallel distributed memory supercomputers.

1. Introduction

Since the advent of modern computers and the first numerical methods for the approximate solution of nonlinear partial differential equations (PDE), the field of computational fluid dynamics (CFD) has been a major driving force for research and new developments in applied mathematics and scientific computing in the last decades, see (e.g., [1,2,3]) for the first finite difference schemes for CFD. Further groundbreaking inventions for the numerical solution of PDE are the finite element (FE) method, which makes use of discrete variational calculus, see (e.g., [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]), and finite volume (FV) schemes [19,20,21,22,23,24,25], which are based on the integral form of the governing PDE system. While finite element schemes are traditionally preferred for computational solid mechanics and stationary elliptic problems, finite volume schemes are widespread in the context of time-dependent hyperbolic partial differential equations, especially when flows become supersonic/supercritical and discontinuous solutions (shock waves) develop, see [26,27].

In [28,29,30,31], a novel family of semi-implicit hybrid FV/FE schemes has been developed and applied to incompressible, weakly compressible, and all Mach number flows, including shock waves. The main idea behind the hybrid FV/FE framework is a flux-vector splitting of the governing PDE system, see [32], where the purely convective (hyperbolic) terms are discretized with an explicit finite volume scheme, while the remaining pressure terms in the flux vector can be rewritten as a second order wave equation and are thus more conveniently discretized at the aid of an implicit finite element method. For the incompressible Navier-Stokes equations, the resulting pressure equation is indeed elliptic, hence the FE approach is clearly more convenient compared to a finite volume discretization. In order to avoid odd-even decoupling and checkerboard effects, the hybrid FV/FE method makes use of a staggered unstructured mesh, where the pressure field is defined in the vertices of a primal simplex mesh composed of triangles or tetrahedra, while all the other flow variables, such as density and momentum, are defined on a dual edge-based/face-based staggered dual mesh, see also [33,34,35]. The use of staggered meshes in CFD goes back to the pioneering work of Harlow and Welch [36] for the solution of the incompressible Navier-Stokes equations and was subsequently also employed, for example, in [37,38]. For a non-exhaustive overview of efficient semi-implicit schemes for CFD, see (e.g., [39,40,41,42,43,44,45,46,47,48,49,50,51,52,53]) and references therein.

In previous work on hybrid FV/FE schemes for the Navier-Stokes equations [28,29,30,31], all simulations have been performed with a classical sequential implementation, allowing the use of only one CPU core for a simulation, even on modern multi-core workstations. In order to reduce the computational cost and perform large simulations with complex geometries and fine meshes needed in CFD for aircraft and vehicles, see (e.g., [54,55,56,57,58]) and references therein for an overview, it is essential to use a parallel implementation. In this paper, the hybrid FV/FE scheme has been parallelized by using the Message-Passing Interface (MPI) standard. The computational domain is partitioned in space via the free domain decomposition software METIS [59], and each processor handles its own data and its own part of the computational domain, hence using a distributed data and distributed workload paradigm. While the parallelization of the explicit finite volume part is straightforward, the implicit discretization of the pressure terms leads to a very large and sparse linear algebraic system that needs to be solved. However, thanks to the underlying basic properties of the governing PDE and its finite element discretization, the resulting discrete pressure system is always symmetric and positive definite. An efficient parallel solution is therefore achieved at the aid of a matrix-free conjugate gradient (CG) method [60]. For more advanced and optimized linear solvers for finite elements, the reader is referred to [61,62,63,64]. The final structure of the matrix-vector product that needs to be evaluated in each CG iteration is similar to the one of an explicit method and can therefore be parallelized as usual at the aid of spatial domain decomposition. In this way, the computational cost is shared among all the processors, both for the explicit FV part as well as for the implicit FE part of the hybrid FV/FE scheme. The parallelism is achieved in all cases via the communication of data among neighboring MPI ranks via shadow/halo zones. For the efficient parallel implementation of explicit schemes for the solution of hyperbolic PDE, the reader is referred to [65,66,67,68,69,70], for example.

By using this parallel implementation of the hybrid FV/FE scheme, numerical simulations involving meshes with several hundreds of millions of elements can be run on modern massively parallel distributed memory supercomputers. Moreover, good parallel performance is obtained, as confirmed by a strong scaling analysis. The numerical solutions obtained with the parallel hybrid FV/FE scheme have been carefully compared with exact analytical, numerical and experimental reference solutions, obtaining in all cases a very good agreement. Currently, the algorithm presented in this paper is limited to simplex meshes, i.e., triangles in 2D and tetrahedra in 3D. Future work will concern the extension to mixed element unstructured meshes, including also prisms, pyramids, and hexahedra, as well as more general polygonal and polyhedral meshes, using e.g., the virtual element method (VEM), see [71,72,73] and references therein. Another limitation of the algorithm is the use of a pure MPI parallelization instead of a mixed MPI + OpenMP implementation.

The structure of this paper is as follows: The mathematical models and the underlying semi-discretization in time that is necessary for the semi-implicit hybrid FV/FE scheme are introduced in Section 2. In particular, the governing PDE systems solved in this paper are the incompressible, the weakly compressible, and the fully compressible Navier-Stokes equations, as well as the shallow water equations for geophysical free surface flows. The spatial discretization of the hybrid FV/FE scheme is presented in Section 3. Section 4 is devoted to describing the parallelization of the staggered semi-implicit hybrid FV/FE approach using the MPI standard. The performance of the parallel scheme is demonstrated via a quantitative speedup analysis. An extensive set of numerical results is presented in Section 5, ranging from the incompressible over the weakly compressible and fully compressible Navier-Stokes equations to the shallow water equations. The paper closes with some concluding remarks and an outlook to future work in Section 6.

2. Mathematical Models and Semi-Discretization in Time

In this section, we recall the mathematical models considered in this paper: the incompressible Navier-Stokes equations, the weakly compressible formulation of the Navier-Stokes equations in terms of the pressure, the fully compressible Navier-Stokes equations in terms of the total energy density, which are suitable for the description of all Mach number flows, and, last but not least, the shallow water equations for the modeling of geophysical free surface flows. Following the hybrid methodology proposed in [28,29,30,31,74], we rely on a flux vector splitting [32] of the associated semi-discrete systems in time. Consequently, we get a first set of transport-diffusion equations that can be efficiently discretized in space using classical explicit finite volume methods. The remaining set of equations can be expressed as a Poisson-type system in the case of the incompressible Navier-Stokes equations and as a second order wave equation for the pressure in the case of the compressible Navier-Stokes equations and the shallow water system. Continuous finite elements are used to solve this system due to the excellent mathematical properties they present for the solution of second order PDE in general and of elliptic problems in particular. Therefore, the corresponding weak formulation of the implicit pressure terms is presented for each mathematical model. Finally, we will recast the resulting systems into a general formulation easing the introduction of the spatial discretization of the equations in the next section.

2.1. Incompressible Navier-Stokes Equations

The incompressible Navier-Stokes equations read

$$\begin{array}{c}\nabla ·{\mathbf{w}}_{\mathbf{v}}=0,\end{array}$$

(1)

$$\begin{array}{c}\frac{\partial {\mathbf{w}}_{\mathbf{v}}}{\partial t}+\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{w}}_{\mathbf{v}}\right)+\nabla p-\nabla ·{\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla \mathbf{v}\right)=\mathbf{0},\end{array}$$

(2)

where ${\mathbf{w}}_{\mathbf{v}}:=\rho \mathbf{v}$ is the linear momentum, $\rho$ the density, $\mathbf{v}$ the velocity vector, ${\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}=\frac{1}{\rho }{\mathbf{w}}_{\mathbf{v}}\otimes {\mathbf{w}}_{\mathbf{v}}$ is the physical flux tensor related to the nonlinear convective terms in the momentum conservation equation, p denotes the pressure, and ${\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}$ corresponds to the viscous stress tensor, which for laminar incompressible flows reads ${\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}=\mu \left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)$, with $\mu$ the dynamic viscosity.

Applying a projection method [53,75], to decouple the computation of the linear momentum unknowns and the pressure field, together with a semi-discretization in time of the above system yields

$$\begin{array}{cc}\hfill {\tilde{\mathbf{W}}}_{\mathbf{v}}& ={\mathbf{W}}_{\mathbf{v}}^n-\Delta t\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{W}}_{\mathbf{v}}^n\right)-\Delta t\nabla {\mathrm{\Pi }}^n+\Delta t\nabla ·{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathbf{W}}_{\mathbf{v}}^{n+1}& ={\tilde{\mathbf{W}}}_{\mathbf{v}}-\Delta t\nabla \left({\mathrm{\Pi }}^{n+1}-{\mathrm{\Pi }}^n\right),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \nabla ·{\mathbf{W}}_{\mathbf{v}}^{n+1}& =0.\hfill \end{array}$$

(5)

Throughout the paper, capital letters are used to identify discrete variables, i.e., ${\mathbf{W}}_{\mathbf{v}}^n$ is the discrete approximation of ${\mathbf{w}}_{\mathbf{v}}\left(\mathbf{x},t^n\right)$, while ${\mathrm{\Pi }}^n=p(\mathbf{x},t^n)$ and so on. Intermediate variables are denoted with a tilde, e.g., ${\tilde{\mathbf{W}}}_{\mathbf{v}}$ is the intermediate variable for the linear momentum including all contributions needed to compute the solution at the new time step but for the spatial gradient of the pressure difference between the two consecutive times, $\nabla \left({\mathrm{\Pi }}^{n+1}-{\mathrm{\Pi }}^n\right)$. Besides, we denote with $\mathrm{\Omega }$ the computational domain and by $\mathrm{\Gamma }:=\partial \mathrm{\Omega }$ its boundary, $\eta$ the outward-pointing normal vector, and by $z\in V_0:=\left\{z\in H^1\left(\mathrm{\Omega }\right):{\int }_{\mathrm{\Omega }}zdV=0\right\}$ a test function. Then, Equations (4) and (5) can be seen as a Poisson-type problem with the following weak formulation:

Weak problem.

Find $\delta {\mathrm{\Pi }}^{n+1}:={\mathrm{\Pi }}^{n+1}-{\mathrm{\Pi }}^n\in V_0$ verifying

$$\begin{array}{c}{\int }_{\mathrm{\Omega }}\nabla \delta {\mathrm{\Pi }}^{n+1}·\nabla z\mathrm{dV}=\frac{1}{\Delta t}{\int }_{\mathrm{\Omega }}{\tilde{\mathbf{W}}}_{\mathbf{v}}^{n+1}·\nabla z\mathrm{dV}-\frac{1}{\Delta t}{\int }_{\mathrm{\Gamma }}{\mathbf{W}}_{\mathbf{v}}^{n+1}·\mathit{\eta }z\mathrm{dS},\forall z\in V_0.\end{array}$$

(6)

2.2. Weakly Compressible Flows

To develop a solver for weakly compressible flows, we consider the mass and momentum equations of the compressible Navier-Stokes system, combined with a pressure equation written in non-conservative form, see [30],

$$\begin{array}{c}\frac{\partial \rho }{\partial t}+\nabla ·{\mathcal{F}}^{\rho }\left(\mathbf{w}\right)=0,\end{array}$$

(7)

$$\begin{array}{c}\frac{\partial {\mathbf{w}}_{\mathbf{v}}}{\partial t}+\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\mathbf{w}\right)+\nabla p-\nabla ·{\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla \mathbf{v}\right)=\rho \mathbf{g},\end{array}$$

(8)

$$\begin{array}{c}\frac{\partial p}{\partial t}+\mathbf{v}·\nabla p-c^2\mathbf{v}·\nabla \rho +c^2\nabla ·\left(\rho \mathbf{v}\right)+\left(\gamma -1\right)\left(\nabla ·\mathbf{q}-{\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla \mathbf{v}\right):\nabla \mathbf{v}\right)=0,\end{array}$$

(9)

with the state vector $\mathbf{w}={(\rho ,\rho \mathbf{v},p)}^T$ and where ${\mathcal{F}}^{\rho }:={\mathbf{w}}_{\mathbf{v}}=\rho \mathbf{v}$ is the flux related to the mass conservation equation, $\mathbf{g}$ is the gravity vector, $\mathbf{q}=-\lambda \nabla \theta$ is the heat flux, $\lambda$ the thermal conductivity coefficient, $\theta$ the temperature, and $c^2={(\partial p/\partial \rho )}_s$ is the square of the isentropic sound speed c. Let us also note that the stress tensor for the laminar compressible regime reads ${\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}=\mu \left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)-\frac{2}{3}\mu \left(\nabla ·\mathbf{v}\right)\mathbf{I}$, with $\mathbf{I}$ the identity matrix, and for the ideal gas equation of state (EOS), we have $c=\sqrt{\frac{\gamma p}{\rho }}$, $p=\rho R\theta$, with R the specific gas constant and $\gamma$ the adiabatic index or ratio of specific heats.

Semi-discretization in time of the former system yields

$$\begin{array}{cc}\hfill {\rho }^{n+1}& ={\rho }^n-\Delta t\nabla ·{\mathcal{F}}^{\rho }\left({\mathbf{W}}^n\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{W}}}_{\mathbf{v}}& ={\mathbf{W}}_{\mathbf{v}}^n-\Delta t\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{W}}^n\right)+\Delta t\nabla ·{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)+\Delta t{\rho }^n\mathbf{g},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\tilde{\mathrm{\Pi }}}_p& =p^n-\Delta t{\mathbf{V}}^n·\nabla {\mathrm{\Pi }}^n-\Delta t\left(\gamma -1\right)\nabla ·{\mathbf{q}}^n,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\tilde{\mathrm{\Pi }}}_{\rho }& =\Delta t{c}^2{\mathbf{V}}^n·\nabla {\rho }^n,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathrm{\Pi }}^{n+1}& =\tilde{\mathrm{\Pi }}-\Delta t{c}^2\nabla ·{\mathbf{W}}_{\mathbf{v}}^{n+1}+\Delta t\left(\gamma -1\right){\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)·\nabla {\mathbf{V}}^n,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\mathbf{W}}_{\mathbf{v}}^{n+1}& ={\tilde{\mathbf{W}}}_{\mathbf{v}}-\Delta t\nabla {\mathrm{\Pi }}^{n+1}\hfill \end{array}$$

(15)

with $\tilde{\mathrm{\Pi }}:={\tilde{\mathrm{\Pi }}}_{\rho }+{\tilde{\mathrm{\Pi }}}_p$. The corresponding second order wave equation for the pressure, derived from (14) and (15), results in

Weak problem.

Find ${\mathrm{\Pi }}^{n+1}\in V_0$ verifying

$$\begin{array}{cc}\hfill \Delta t{}^2{\int }_{\mathrm{\Omega }}\nabla {\mathrm{\Pi }}^{n+1}·\nabla z\mathrm{dV}+\frac{1}{c^2}{\int }_{\mathrm{\Omega }}{\mathrm{\Pi }}^{n+1}z\mathrm{dV}=& \Delta t{\int }_{\mathrm{\Omega }}{\tilde{\mathbf{W}}}_{\mathbf{v}}·\nabla z\mathrm{dV}-\Delta t{\int }_{\mathrm{\Gamma }}{\mathbf{W}}_{\mathbf{v}}^{n+1}·\eta z\mathrm{dS}\hfill \\ \hfill & +\frac{1}{c^2}{\int }_{\mathrm{\Omega }}\tilde{\mathrm{\Pi }}z\mathrm{dV}-\frac{\Delta t}{c^2}{\int }_{\mathrm{\Omega }}\left(\gamma -1\right){\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)·\nabla {\mathbf{V}}^{n}z\mathrm{dV},\forall z\in V_0.\hfill \end{array}$$

(16)

2.3. Compressible Navier-Stokes Equations for All Mach Number Flows

To get an all Mach number flow solver, we go back to the classical conservative formulation of the compressible Navier-Stokes equations in terms of the total energy density, $w_E:=\rho E$ [31],

$$\begin{array}{c}\frac{\partial \rho }{\partial t}+\nabla ·{\mathcal{F}}^{\rho }\left(\mathbf{w}\right)=0,\end{array}$$

(17)

$$\begin{array}{c}\frac{\partial {\mathbf{w}}_{\mathbf{v}}}{\partial t}+\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\mathbf{w}\right)+\nabla p-\nabla ·{\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla \mathbf{v}\right)=\rho \mathbf{g},\end{array}$$

(18)

$$\begin{array}{c}\frac{\partial w_E}{\partial t}+\nabla ·\left[\frac{1}{\rho }{\mathbf{w}}_{\mathbf{v}}\left(w_E+p\right)\right]-\nabla ·\left(\frac{1}{\rho }{\mathit{\tau }}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla \mathbf{v}\right)·{\mathbf{w}}_{\mathbf{v}}\right)+\nabla ·\mathbf{q}=\mathbf{g}·{\mathbf{w}}_{\mathbf{v}}.\end{array}$$

(19)

Here, the state vector is $\mathbf{w}=(\rho ,\rho v,\rho E)$.

To discretize the compressible system, we follow the flux splitting techniques introduced in [32,48], obtaining the following semi-discrete system in time:

$$\begin{array}{cc}\hfill W_{\rho }^{n+1}& =W_{\rho }^n-\Delta t\nabla ·{\mathcal{F}}^{\rho }\left({\mathbf{W}}^n\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{W}}}_{\mathbf{v}}& ={\mathbf{W}}_{\mathbf{v}}^n-\Delta t\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{W}}^n\right)-\Delta t\nabla {\mathrm{\Pi }}^n+\Delta t\nabla ·{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)+\Delta t{\rho }^n\mathbf{g},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\tilde{W}}_E& =W_E^n-\Delta t\nabla {\mathcal{F}}^K\left({\mathbf{W}}^n\right)+\Delta t\nabla ·\left(\frac{1}{{\rho }^n}{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right){\mathbf{W}}_{\mathbf{v}}^n\right)+\Delta t\nabla ·{\mathbf{Q}}^n+\Delta t\mathbf{g}·{\mathbf{W}}_{\mathbf{v}}^n,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{1}{\gamma -1}\delta {\mathrm{\Pi }}^{n+1,k+1}& -\Delta t{}^2\nabla ·\left(H^{n+1,k}\nabla \delta {\mathrm{\Pi }}^{n+1,k+1}\right)={\tilde{W}}_E-\frac{{\mathrm{\Pi }}^n}{\gamma -1}-{\left(\rho K\right)}^{n+1,k}-\Delta t\nabla ·\left(H^{n+1,k}{\tilde{\mathbf{W}}}_{\mathbf{v}}\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathbf{W}}_{\mathbf{v}}^{n+1,k+1}& ={\tilde{\mathbf{W}}}_{\mathbf{v}}-\Delta t\nabla \delta {\mathrm{\Pi }}^{n+1,k+1},\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathrm{\Pi }}^{n+1,k+1}& ={\mathrm{\Pi }}^n+\delta {\mathrm{\Pi }}^{n+1,k+1},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill W_E^{n+1}& ={\tilde{W}}_E-\Delta t\nabla ·\left(H^{n+1,k+1}{\mathbf{W}}_{\mathbf{v}}^{n+1,k+1}\right),\hfill \end{array}$$

(26)

where ${\mathcal{F}}^K:=K^n{\mathbf{W}}_{\mathbf{v}}^n$ is the flux related to the kinetic energy. Let us remark that in (23)–(26), we have introduced a Picard iteration procedure, k the Picard index, in order to avoid the solution of a highly nonlinear pressure system, further details can be found in [31]. The weak formulation related to the pressure subsystem (23) reads

Weak problem.

Find $\delta {\mathrm{\Pi }}^{n+1,k+1}={\mathrm{\Pi }}^{n+1,k+1}-{\mathrm{\Pi }}^n\in V_0$ verifying

$$\begin{array}{cc}\hfill \frac{1}{\gamma -1}{\int }_{\mathrm{\Omega }}\delta {\mathrm{\Pi }}^{n+1,k+1}z\mathrm{dV}& -\Delta t{}^2{\int }_{\mathrm{\Omega }}{H}^{n+1,k}\nabla \delta {\mathrm{\Pi }}^{n+1,k+1}·\nabla z\mathrm{dV}=-\Delta t{\int }_{\mathrm{\Gamma }}{H}^{n+1,k}{\mathbf{W}}_{\mathbf{v}}^{n+1}·\eta z\mathrm{dS}\hfill \\ \hfill & +{\int }_{\mathrm{\Omega }}\left({\tilde{W}}_E-{\left(\rho K\right)}^{n+1,k}\right)z\mathrm{dV}-\frac{1}{\gamma -1}{\int }_{\mathrm{\Omega }}{\mathrm{\Pi }}^{n}z\mathrm{dV}+\Delta t{\int }_{\mathrm{\Omega }}{H}^{n+1,k}{\tilde{\mathbf{W}}}_{\mathbf{v}}·\nabla z\mathrm{dV},\forall z\in V_0.\hfill \end{array}$$

(27)

2.4. Shallow Water Equations

Following [40], the well-known shallow water system of partial differential equations can be written as

$$\begin{array}{ccc}\hfill \frac{\partial \eta }{\partial t}+\nabla ·{\mathbf{w}}_{\mathbf{v}}& =& 0,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \frac{\partial {\mathbf{w}}_{\mathbf{v}}}{\partial t}+\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\mathbf{w}\right)+gh\nabla \eta & =& -\gamma {\mathbf{w}}_{\mathbf{v}},\hfill \end{array}$$

(29)

where $\eta :=h+b$ is the free surface elevation, h is the water depth, b is the bottom elevation, that we assume constant in time, ${\mathbf{w}}_{\mathbf{v}}=h\mathbf{v}$ is the conservative variable related to the velocity field, ${\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}:=\frac{1}{h}{\mathbf{w}}_{\mathbf{v}}\otimes {\mathbf{w}}_{\mathbf{v}}$ is the flux of the momentum conservation equation, g is the gravity constant, and $\gamma$ is a shorthand notation for the bottom friction.

$$\gamma =\frac{g}{{\kappa }_s^2{h}^{\frac{4}{3}}}\left|\mathbf{v}\right|,$$

(30)

with ${\kappa }_s=n^{-1}$ the Strickler coefficient and n the Manning coefficient. An important aspect that should be noticed is that the hybrid FV/FE algorithm proposed is well balanced by construction so water at rest solutions are preserved [33,76,77].

After discretization in time and some manipulations of the resulting equations, see [74] for further details, we get

$$\begin{array}{cc}\hfill {\tilde{\mathbf{W}}}_{\mathbf{v}}& ={\mathbf{W}}_{\mathbf{v}}^n-\Delta t\nabla ·{\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{W}}^n\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\eta }^{n+1}-\Delta t{}^{2}g\nabla ·\left(\frac{H^n}{1+\Delta t{\gamma }^n}\nabla {\eta }^{n+1}\right)& ={\eta }^n-\Delta t\nabla ·\left(\frac{1}{1+\Delta t{\gamma }^n}{\tilde{\mathbf{W}}}_{\mathbf{v}}\right),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathbf{W}}_{\mathbf{v}}^{n+1}& =\frac{1}{1+\Delta t{\gamma }^n}\left({\tilde{\mathbf{W}}}_{\mathbf{v}}-\Delta tg{H}^n\nabla {\eta }^{n+1}\right).\hfill \end{array}$$

(33)

Similar to what we have already seen for the pressure subsystems for the compressible Navier-Stokes equations, we get a second order wave equation for the unknown free surface. Besides, to increase the accuracy in time, we employ the $\theta$-method that reduces to the classical Crank-Nicolson scheme for $\theta =0.5$. The final weak formulation results in

Weak problem.

Find ${\eta }^{n+1}\in V_0$ verifying

$$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{\eta }^{n+1}z\mathrm{dV}& +\Delta t{}^2{\theta }^{2}g{\int }_{\mathrm{\Omega }}\left(\frac{H^n}{1+\Delta t{\gamma }^n}\nabla {\eta }^{n+1}\right)·\nabla z\mathrm{dV}={\int }_{\mathrm{\Omega }}{\eta }^{n}zdV+\Delta t{\int }_{\mathrm{\Omega }}{\widetilde{\tilde{\mathbf{W}}}}_{\mathbf{v}}\nabla z\mathrm{dV}\hfill \\ \hfill & -\Delta t{\int }_{\mathrm{\Omega }}{\mathbf{W}}_{\mathbf{v}}^{n+\theta }·\mathbf{n}z\mathrm{dV}-\Delta t{}^2\theta \left(1-\theta \right)g{\int }_{\mathrm{\Omega }}\frac{H^n}{1+\Delta t{\gamma }^n}\nabla {\eta }^n·\nabla z\mathrm{dV},\forall z\in V_0.\hfill \end{array}$$

(34)

In (34), we have introduced a new intermediate variable

$${\widetilde{\tilde{\mathbf{W}}}}_{\mathbf{v}}=\left(1-\theta \right){\mathbf{W}}_{\mathbf{v}}^n+\frac{\theta }{1+\Delta t{\gamma }^n}{\tilde{\mathbf{W}}}_{\mathbf{v}}.$$

2.5. General Formulation

As mentioned before, the former systems can be recast into a general formulation which eases the development of a general code for its solution. In particular, we consider two different sets of equations: the first one is related to the transport-diffusion equations, while the second one corresponds to the second order pressure system. Moreover, an extra set of equations allows the final update of the intermediate variables computed from the transport-diffusion equations once the implicit pressure system has been solved.

Denoting $\mathbf{Q}$ the vector of conservative variables, the general form of the transport-diffusion subsystems reads:

$$\tilde{\mathbf{Q}}={\mathbf{Q}}^n-\Delta t\nabla ·\mathcal{F}\left({\mathbf{Q}}^n\right)-\Delta t\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n-\Delta t\nabla {\mathbf{P}}^n+\Delta t\nabla ·\mathbf{T}\left(\nabla {\mathbf{Q}}^n\right)+\Delta t{\mathbf{f}}_1\left(\nabla {\mathbf{Q}}^n\right)+\Delta t{\mathbf{f}}_2\left({\mathbf{Q}}^n\right),$$

(35)

while the second order pressure equation can be reshaped into:

$$\begin{array}{c}{\int }_{\mathrm{\Omega }}{G}_1\left({\mathbf{Q}}^n\right)P^{n+1}z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{H}_1\left({\mathbf{Q}}^{n+1,k},{\mathbf{Q}}^n\right)\nabla P^{n+1}·\nabla z\mathrm{dV}=-{\int }_{\mathrm{\Gamma }}{\mathbf{H}}_2\left({\mathbf{Q}}^{n+1,k},{\mathbf{Q}}^{n+1},{\mathbf{Q}}^{n+\theta }\right)·\mathbf{n}z\mathrm{dS}\\ +{\int }_{\mathrm{\Omega }}{G}_2\left(P^n,\tilde{P}\right)zdV+{\int }_{\mathrm{\Omega }}{\mathbf{H}}_3\left({\mathbf{Q}}^{n+1,k},\tilde{\mathbf{Q}}\right)·\nabla z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{G}_3\left({\mathbf{Q}}^{n+1,k},\tilde{\mathbf{Q}},{\mathbf{Q}}^n\right)z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{\mathbf{G}}_4\left({\mathbf{Q}}^n,\nabla P^n\right)·\nabla z\mathrm{dV}.\end{array}$$

(36)

The different terms in the former equation have been clustered into the auxiliary functions $G_1,G_2,G_3,{\mathbf{G}}_4,H_1,{\mathbf{H}}_2,{\mathbf{H}}_3$ to ease the introduction of the numerical method in the next sections. Note that not all arguments are necessary for all models. In what follows, we detail the values taken by each term in (35)–(36) according to the four models considered:

• Incompressible Navier-Stokes equations:

  $$\begin{array}{c}\tilde{\mathbf{Q}}={\tilde{\mathbf{W}}}_{\mathbf{v}},{\mathbf{Q}}^n={\mathbf{W}}_{\mathbf{v}}^n,\mathcal{F}\left({\mathbf{Q}}^n\right)={\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{Q}}^n\right),{\mathbf{P}}^n={\mathrm{\Pi }}^n,\\ \mathbf{T}\left(\nabla {\mathbf{Q}}^n\right)={\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right),\mathbf{B}\left({\mathbf{Q}}^n\right)={\mathbf{f}}_1\left(\nabla {\mathbf{Q}}^n\right)={\mathbf{f}}_2\left({\mathbf{Q}}^n\right)=\mathbf{0},\\ P^{n+1,k+1}=\delta {\mathrm{\Pi }}^{n+1},G_1=0,H_1=1,{\mathbf{H}}_2=\frac{1}{\Delta t}{\mathbf{W}}_{\mathbf{v}}^{n+1},{\mathbf{H}}_3=\frac{1}{\Delta t}{\tilde{\mathbf{W}}}_{\mathbf{v}},G_2=G_3=0,{\mathbf{G}}_4=\mathbf{0}.\end{array}$$
• Weakly compressible flows:

  $$\begin{array}{c}\tilde{\mathbf{Q}}=\left(\begin{array}{c}{\rho }^{n+1}\\ {\tilde{\mathbf{W}}}_{\mathbf{v}}\\ {\tilde{\mathrm{\Pi }}}_p\end{array}\right),{\mathbf{Q}}^n=\left(\begin{array}{c}{\rho }^n\\ {\mathbf{W}}_{\mathbf{v}}^n\\ {\mathrm{\Pi }}^n\end{array}\right),\mathcal{F}\left({\mathbf{Q}}^n\right)=\left(\begin{array}{c}{\mathcal{F}}^{\rho }\left({\mathbf{Q}}^n\right)\\ {\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{Q}}^n\right)\\ 0\end{array}\right),\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n=\left(\begin{array}{c}0\\ \mathbf{0}\\ {\mathbf{V}}^n·\nabla {\mathrm{\Pi }}^n\end{array}\right),{\mathbf{P}}^n=\mathbf{0},\\ \mathbf{T}\left(\nabla {\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ {\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)\\ 0\end{array}\right),{\mathbf{f}}_1\left(\nabla {\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ \mathbf{0}\\ \left(\gamma -1\right)\nabla ·{\mathbf{Q}}^n\end{array}\right),{\mathbf{f}}_2\left({\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ {\rho }^n\mathbf{g}\\ 0\end{array}\right),\tilde{\mathrm{\Pi }}={\tilde{\mathrm{\Pi }}}_p+\Delta t{c}^2{\mathbf{v}}^n·\nabla {\rho }^n,\\ P^{n+1,k+1}={\mathrm{\Pi }}^{n+1},G_1=\frac{1}{c^2},H_1=\Delta t{}^2,{\mathbf{H}}_2=\Delta t{\mathbf{W}}_{\mathbf{v}}^{n+1},{\mathbf{H}}_3=\Delta t{\tilde{\mathbf{W}}}_{\mathbf{v}},\\ G_2=\frac{1}{c^2}\tilde{\mathrm{\Pi }},G_3=\frac{\Delta t}{c^2}{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)·\nabla {\mathbf{V}}^n,{\mathbf{G}}_4=\mathbf{0}.\end{array}$$
  Let us remark that in this case, an additional term, ${\mathrm{\Pi }}_{\rho }$, is computed between the approximation of the transport equations and the second order wave equation.
• Compressible Navier-Stokes equations for all Mach number flows:

  $$\begin{array}{c}\tilde{\mathbf{Q}}=\left(\begin{array}{c}{\rho }^{n+1}\\ {\tilde{\mathbf{W}}}_{\mathbf{v}}\\ {\tilde{W}}_E\end{array}\right),{\mathbf{Q}}^n=\left(\begin{array}{c}{\rho }^n\\ {\mathbf{W}}_{\mathbf{v}}^n\\ W_E^n\end{array}\right),\mathcal{F}\left({\mathbf{Q}}^n\right)=\left(\begin{array}{c}{\mathcal{F}}^{\rho }\left({\mathbf{Q}}^n\right)\\ {\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{Q}}^n\right)\\ {\mathcal{F}}^K\left({\mathbf{Q}}^n\right)\end{array}\right),\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n=\mathbf{0},\\ \nabla {\mathbf{P}}^n=\left(\begin{array}{c}0\\ \nabla {\mathrm{\Pi }}^n\\ 0\end{array}\right),\mathbf{T}\left(\nabla {\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ {\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right)\\ \frac{1}{{\rho }^n}{\mathbf{T}}^{{\mathbf{w}}_{\mathbf{v}}}\left(\nabla {\mathbf{V}}^n\right){\mathbf{W}}_{\mathbf{v}}^n\end{array}\right),{\mathbf{f}}_1\left(\nabla {\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ \mathbf{0}\\ \nabla ·{\mathbf{Q}}^n\end{array}\right),{\mathbf{f}}_2\left({\mathbf{Q}}^n\right)=\left(\begin{array}{c}0\\ {\rho }^n\mathbf{g}\\ \mathbf{g}·{\mathbf{W}}_{\mathbf{v}}^n\end{array}\right),\\ P^{n+1,k+1}=\delta {\mathrm{\Pi }}^{n+1,k+1},G_1=\frac{1}{\gamma -1},H_1=-\Delta t{}^2{H}^{n+1,k},{\mathbf{H}}_2=\Delta t{H}^{n+1,k}{\mathbf{W}}_{\mathbf{v}}^{n+1},\\ {\mathbf{H}}_3=\Delta t{H}^{n+1,k}{\tilde{\mathbf{W}}}_{\mathbf{v}},G_2=-\frac{1}{\gamma -1}{\mathrm{\Pi }}^n,G_3={\tilde{W}}_E-{\left(\rho K\right)}^{n+1,k},{\mathbf{G}}_4=\mathbf{0}.\end{array}$$
• Shallow water equations:

  $$\begin{array}{c}\tilde{\mathbf{Q}}={\tilde{\mathbf{W}}}_{\mathbf{v}},{\mathbf{Q}}^n={\mathbf{W}}_{\mathbf{v}}^n,\mathcal{F}\left({\mathbf{Q}}^n\right)={\mathcal{F}}^{{\mathbf{w}}_{\mathbf{v}}}\left({\mathbf{Q}}^n\right),{\mathbf{P}}^n=\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n=\mathbf{T}\left({\mathbf{Q}}^n\right)={\mathbf{f}}_2\left(\nabla {\mathbf{Q}}^{n+1}\right)={\mathbf{f}}_2\left({\mathbf{Q}}^{n+1}\right)=\mathbf{0},\\ P^{n+1,k+1}={\eta }^{n+1},G_1=1,H_1=\Delta t{}^2{\theta }^{2}g\frac{H^n}{1+\Delta t{\gamma }^n},{\mathbf{H}}_2=\Delta t{\mathbf{W}}_{\mathbf{v}}^{n+\theta },{\mathbf{H}}_3=\Delta t{\widetilde{\tilde{\mathbf{W}}}}_{\mathbf{v}},\\ G_2={\eta }^n,G_3=0,{\mathbf{G}}_4=-\Delta t{}^2\theta \left(1-\theta \right)g\frac{H^n}{1+\Delta t{\gamma }^n}\nabla {\eta }^n.\end{array}$$

3. Spatial Discretization

To complete the description of the numerical method, we now present the spatial discretization of the system (35)–(36). The hybrid FV/FE algorithms [28,29,30,31,74] can be organized into four steps:

• Transport-diffusion stage. A first intermediate approximation is given for the conserved variables by solving (35) using explicit finite volume methods.
• Interpolation stage. For spatial discretization, two staggered grids are employed for each numerical method, so interpolation of some variables between the two different grids is needed.
• Second order pressure equation stage. System (36) is approximated using implicit continuous finite elements.
• Correction stage. The intermediate approximation of the conservative variables obtained in the first stage is updated with the solution of the second order pressure system.

In what follows, we first introduce the staggered unstructured grids employed in our hybrid FV/FE framework. Then, we focus on the finite volume method applied to the transport-diffusion equations. Finally, some key points of the finite element methodology and the correction stage are also provided.

3.1. Staggered Unstructured Grid

In order to deal with complex geometries, unstructured simplex meshes are used to discretize the computational domain $\mathrm{\Omega }$. In particular, edge-type (2D)/face-type (3D) unstructured staggered meshes are considered in the spatial discretization. To introduce the notation needed to explain the numerical algorithm, we summarize here the mesh structure in 2D, further details on the generation of three-dimensional face-based staggered meshes can be found in [28,47,50].

The spatial computational domain $\mathrm{\Omega }\in {\mathbb{R}}^2$ is covered by a set of disjoint triangular elements $\left\{T_k,k=1,\dots ,n_{el}\right\}$, with vertex $V_j$, $j=\left\{1,\dots ,n_{ver}\right\}$, which form the primal mesh, see the left sketch in Figure 1. Moreover, we denote by B and $B^{\prime }$ the barycentres of two primal elements sharing one interior edge. To build the dual finite volume mesh, two different kinds of elements are defined: the interior dual elements, which are built by merging the two triangles which result from connecting this common edge with their barycentres (see Figure 1 centre), and the dual boundary elements which are built connecting the two vertices of a boundary edge in the primal mesh to the two barycentres of the triangles sharing this common edge. Each dual element $C_i$ is also called cell. On this dual mesh, the nodes $\left\{N_i,i=1,\dots ,n_{nod}\right\}$ are the barycentres of the edges of the elements on the primal mesh. Similarly, in 3D the primal elements are tetrahedra. Each interior dual element is built by merging the two tetrahedra created connecting one interior face with the barycentres of the two primal elements which contain this face, and each boundary dual element is defined by the vertices of one boundary face and the barycentre of the tetrahedron in the primal mesh containing that face. Thus, the nodes $N_i$ are the barycentres of the faces of the primal elements. Further details on the 2D and 3D staggered meshes are shown in Figure 1 and Figure 2, respectively. The remaining mesh-related notation to be used is as follows:

• ${\mathcal{K}}_i$ is the set of neighboring nodes of a node $N_i$, which are the barycentres of the dual cells sharing an edge/face with $C_i$.
• ${\mathrm{\Gamma }}_i$ is the boundary of cell $C_i$ and ${\tilde{\mathit{\eta }}}_i$ is its outward unit normal.
• ${\mathrm{\Gamma }}_{ij}$ is the edge/face between cells $C_i$, and $C_j$. $N_{ij}$ is the barycentre of ${\mathrm{\Gamma }}_{ij}$ (Figure 1 right), and ${\displaystyle {\mathrm{\Gamma }}_i=\bigcup _{N_j\in {\mathcal{K}}_i}{\mathrm{\Gamma }}_{ij}}$.
• $\left|C_i\right|$ is the area/volume of cell $C_i$.
• ${\tilde{\mathit{\eta }}}_{ij}$ is the outward unit normal vector to ${\mathrm{\Gamma }}_{ij}$. We define $\mathit{\eta }{}_{ij}:={\tilde{\mathit{\eta }}}_{ij}\left|\right|\mathit{\eta }{}_{ij}\left|\right|$, where, $\left|\right|\mathit{\eta }{}_{ij}\left|\right|$ represents the length/area of ${\mathrm{\Gamma }}_{ij}$.

3.2. Transport-Diffusion Stage

Integrating (35) on each dual control volume $C_i$, applying the Gauss theorem to advection, pressure, and diffusion terms, and considering a path conservative approach for the non-conservative products yields

$$\begin{array}{cc}\hfill {\tilde{\mathbf{Q}}}_i=& {\mathbf{Q}}_i^n-\frac{\Delta t}{\left|C_i\right|}{\int }_{{\mathrm{\Gamma }}_i}\mathcal{F}\left({\mathbf{Q}}^n\right)\tilde{\mathit{\eta }}{}_i\mathrm{dS}-\frac{\Delta t}{\left|C_i\right|}{\int }_{{\mathrm{\Gamma }}_i}\mathcal{D}\left({\mathbf{Q}}^n\right)\tilde{\mathit{\eta }}{}_i\mathrm{dS}-\frac{\Delta t}{\left|C_i\right|}{\int }_{C_i}\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n\mathrm{dV}-\frac{\Delta t}{\left|C_i\right|}{\int }_{{\mathrm{\Gamma }}_i}{\mathbf{P}}^n\tilde{\mathit{\eta }}{}_i\mathrm{dS}\hfill \\ \hfill & +\frac{\Delta t}{\left|C_i\right|}{\int }_{{\mathrm{\Gamma }}_i}\mathbf{T}\left(\nabla {\mathbf{Q}}^n\right)\tilde{\mathit{\eta }}{}_i\mathrm{dS}+\frac{\Delta t}{\left|C_i\right|}{\int }_{C_i}{\mathbf{f}}_1\left(\nabla {\mathbf{Q}}^n\right)\mathrm{dV}+\frac{\Delta t}{\left|C_i\right|}{\int }_{C_i}{\mathbf{f}}_2\left({\mathbf{Q}}^n\right)\mathrm{dV}.\hfill \end{array}$$

(37)

To approximate the flux term and get a second order scheme in space and time, we make use of the Rusanov flux function [78], combined with a local ADER methodology [29,79]. We first split the integral on the cell boundary in the sum of the integrals at the cell edges/faces. For an arbitrary conservative variable Q, the corresponding numerical flux function reads

$$\begin{array}{c}{\displaystyle \varphi \left({\overline{Q}}_i^n,{\overline{Q}}_j^n,\mathit{\eta }{}_{ij}\right)=\frac{1}{2}\left({\mathcal{Z}}^Q\left({\overline{\mathbf{Q}}}_i^n,\mathit{\eta }{}_{ij}\right)+{\mathcal{Z}}^Q\left({\overline{\mathbf{Q}}}_j^n,\mathit{\eta }{}_{ij}\right)\right)-\frac{1}{2}{\tilde{\alpha }}_{RS,ij}^n\left({\overline{Q}}_j^n-{\overline{Q}}_i^n\right)}\end{array}$$

(38)

with ${\mathcal{Z}}^Q({\overline{\mathbf{Q}}}^n,\mathit{\eta }{}_i):={\mathcal{F}}^Q\left({\overline{\mathbf{Q}}}^n\right)\mathit{\eta }{}_i$ the physical flux in the normal direction and

$${\tilde{\alpha }}_{RS,ij}^n={\tilde{\alpha }}_{RS}^n\left({\overline{\mathbf{Q}}}_i^n,{\overline{\mathbf{Q}}}_j^n,\mathit{\eta }{}_{ij}\right):={\alpha }_{RS}^n\left({\overline{\mathbf{Q}}}_i^n,{\overline{\mathbf{Q}}}_j^n,\mathit{\eta }{}_{ij}\right)+c_{\alpha }∥\mathit{\eta }{}_{ij}∥$$

(39)

the Rusanov coefficient, where ${\alpha }_{RS}^n$ is the maximum signal speed on the edge, computed as the maximum absolute eigenvalue of the complete transport-diffusion subsystem, and $c_{\alpha }\in {\mathbb{R}}_0^{+}$ is an optional artificial viscosity coefficient that can be activated to improve the stability properties of the scheme. The value of the half in time evolved reconstructed variables ${\overline{\mathbf{Q}}}^n$ can be obtained following the classical ADER procedure:

Step 1.

Piecewise polynomial reconstruction of the data in the neighboring of each edge/face of the cell. To circumvent Godunov’s theorem and obtain a stable scheme, the gradients for the definition of the polynomial are selected using a nonlinear Essentially Non-Oscillatory (ENO) methodology [80]. As alternative limiting strategies, we employ the classical mimmod limiter [23,81] and the more restrictive Barth-Jespersen limiter, see [82]. Besides, the computation of the gradients is done using a Galerkin approach which benefits from the staggered meshes arrangement to get a reduced stencil [30].

Step 2.

Computation of the boundary-extrapolated values in the barycentre of the faces by evaluating the piecewise polynomials.

Step 3.

Half in time evolution of the data by using a temporal Taylor series expansion. The time derivatives are replaced by spatial derivatives considering a Cauchy-Kovalevskaya procedure based on the original model, i.e., on the continuous PDEs where the splitting technique has not been applied yet.

The second and third integrals in (37) correspond to the non-conservative product computed following the path conservative schemes introduced in [76,83,84,85,86]. Consequently, its contribution is divided into a first part accounting for the jumps between cells and the smooth contribution from inside each cell necessary for getting the second order scheme. For the jump term, needed to preserve equilibrium solutions of the PDEs exactly at the discrete level, we perform the integral on the face following [76,85,87]:

$$\frac{\Delta t}{\left|C_i\right|}{\int }_{{\mathrm{\Gamma }}_{ij}}\mathcal{D}\left({\mathbf{Q}}^n\right)\tilde{\mathit{\eta }}{}_{ij}\mathrm{dS}\approx \frac{1}{2}\sum _{j\in {\mathcal{K}}_i}{\tilde{\mathbf{B}}}_{ij}\left({\mathbf{Q}}_j^n-{\mathbf{Q}}_i^n\right),{\tilde{\mathbf{B}}}_{ij}={\int }_0^1\mathbf{B}\left(\psi \left({\mathbf{Q}}_i^n,{\mathbf{Q}}_j^n,s\right)\right)·\tilde{\mathit{\eta }}{}_{ij}\mathrm{ds}.$$

(40)

In particular, $\psi =\psi \left({\mathbf{Q}}_i^n,{\mathbf{Q}}_j^n,s\right)={\mathbf{Q}}_i^n+s\left({\mathbf{Q}}_j^n-{\mathbf{Q}}_i^n\right)$ has been chosen to be the simple straight line segment path, and the trapezoidal rule is used to evaluate the integral in the definition of the Roe-type matrix ${\tilde{\mathbf{B}}}_{ij}$. Meanwhile, the smooth part of the non-conservative product,

$$\frac{\Delta t}{\left|C_i\right|}{\int }_{C_i}\mathbf{B}\left({\mathbf{Q}}^n\right)·\nabla {\mathbf{Q}}^n\mathrm{dV},$$

(41)

is approximated, assuming a constant value by cell. First, the gradients in the dual elements are obtained as a weighted average of the values at the primal elements computed using the Galerkin approach. Then, we perform the product of matrix $\mathbf{B}\left({\mathbf{Q}}_i^n\right)$ by the obtained ${\left(\nabla {\mathbf{Q}}^n\right)}_{|C_i}$.

Regarding the pressure term, we assume a known value at each dual vertex and approximate the integral on the face using the corresponding quadrature rule. In case the vertex corresponds to the barycentre of a primal element, the pressure value is taken as the average value at the vertex of the element.

Viscous terms are computed by considering the normal projection of the Galerkin approximation of the gradients on each face. Again a Taylor series expansion in time jointly with a Cauchy-Kovalevskaya procedure, in which advection terms are neglected, are employed to get a second order scheme.

Finally, algebraic source terms can be computed using a high enough quadrature rule. In case source terms depend on the system unknowns, a constant value per cell is assumed while Gauss theorem is applied when its gradients are involved.

3.3. Interpolation Stage

Since two staggered grids are used, interpolation of data between them becomes necessary. To minimize artificial diffusion due to this process, we assume the averaged quantity on dual cells corresponding with the value on the node, and we pass to primal cells following

$${\mathbf{Q}}_k=\frac{1}{lnf}\sum _{i\in {\mathcal{K}}_k}{\mathbf{Q}}_i$$

(42)

with ${\mathcal{K}}_k$ the set of indexes pointing to the faces of the primal element $T_k$, $C_i$ the dual cell related to face i, and $lnf$ the local number of vertex of the faces of the primal element, namely, $lnf=3$ in 2D and $lnf=4$ in 3D. The proposed interpolation leads to a provable exact interpolation technique to pass data between the two different staggered grids, see [88].

On the other hand, the interpolation back from the primal vertices to the dual cells is achieved as

$${\mathbf{Q}}_i=\frac{1}{lnvf}\sum _{k\in {\mathcal{V}}_i}{\mathbf{Q}}_k,$$

(43)

where ${\mathcal{V}}_i$ contains the index of the vertex on the face containing the node $N_i$ and $lnvf$ is the number of vertex at each face. Further interpolations considered, including the computation of the data in the vertex of primal elements and calculation of auxiliary variables from mixed data originally coming from both grids, can be found, e.g., in [31].

3.4. Finite Element Method

As pointed out in Section 2.5, all second order pressure equations of the considered mathematical models can be recast into a general formulation whose corresponding weak problem reads:

Weak problem.

Find $P^{n+1}\in V_0$ verifying

$$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{G}_1\left({\mathbf{Q}}^n\right)P^{n+1}z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{H}_1\left({\mathbf{Q}}^{n+1,k},{\mathbf{Q}}^n\right)\nabla P^{n+1}·\nabla z\mathrm{dV}=-{\int }_{\mathrm{\Gamma }}{\mathbf{H}}_2\left({\mathbf{Q}}^{n+1,k},{\mathbf{Q}}^{n+1},{\mathbf{Q}}^{n+\theta }\right)·\mathbf{n}z\mathrm{dS}\\ \hfill +{\int }_{\mathrm{\Omega }}{G}_2\left(P^n,\tilde{P}\right)zdV+{\int }_{\mathrm{\Omega }}{\mathbf{H}}_3\left({\mathbf{Q}}^{n+1,k},\tilde{\mathbf{Q}}\right)·\nabla z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{G}_3\left({\mathbf{Q}}^{n+1,k},\tilde{\mathbf{Q}},{\mathbf{Q}}^n\right)z\mathrm{dV}+{\int }_{\mathrm{\Omega }}{\mathbf{G}}_4\left({\mathbf{Q}}^n,\nabla P^n\right)·\nabla z\mathrm{dV}.\end{array}$$

(44)

Discretization in space at the aid of ${\mathbb{P}}^1$ Lagrange finite elements will provide the solution $P^{n+1}$ at the vertices of the primal elements. The corresponding linear algebraic system is symmetric and positive definite and can be efficiently solved using a matrix-free conjugate gradient method [60]. Let us remark that for fully compressible flows, a Picard iteration procedure is applied to avoid the complex nonlinear system resulting from the presence of both the enthalpy and the kinetic energy in the integrals on the left-hand side of (36), see [31] for further details. Besides, the $\theta$-method, introduced in [74] for the shallow water case, could be applied to increase the time accuracy of the overall scheme.

3.5. Correction Stage

Since the intermediate approximations obtained in the transport diffusion stage lack for the contributions of $P^{n+1}$ related terms, a correction step must be performed at the end of each time step. To this end, in the correction stage or post-projection stage, we solve an extra update equation of the form

$${\mathbf{Q}}_i^{n+1}={\mathbf{C}}_1\left({\mathbf{Q}}^n,P^n\right){\tilde{\mathbf{Q}}}_i-\Delta t{\mathbf{C}}_2\left({\mathbf{Q}}^n,P^n,P^{n+1},\nabla P^{n+1}\right)$$

(45)

with the functions ${\mathbf{C}}_1,{\mathbf{C}}_2$ defined according to Equations (4), (15), (24)–(26) and (33), for each mathematical model considered. Note that ${\mathbf{C}}_2$ depends not only on the value of the variables at the new time step but also on their gradients. To compute this contribution, the gradient is first computed at each primal element, and then a weighted average provides the value in the dual elements.

4. Parallel Implementation

In [28,29,30,31,74], the hybrid FV/FE approach has been widely described for incompressible, weakly compressible, all Mach number flows, and shallow water equations showing several results run using a sequential code, written in Fortran language. There exist many circumstances that lead us to have a high computational cost, such as a large computational domain, an extremely fine mesh, a complex geometry, or a large final time, and it, therefore, becomes necessary to implement the numerical scheme in parallel, hence allowing us to reduce the computational cost of the simulation. High-Performance Computing (HPC) involves different techniques to achieve efficient and fast simulations, such as distributed or shared memory parallelisms, vectorization, or memory access optimization. For the family of hybrid FV/FE methods, we consider parallel programming on distributed memory environments with the Message-Passing Interface (MPI) standard [89], to distribute data and workload. The proposed parallel implementation is based on decomposing the computational domain into many subdomains allotted to various processors (MPI ranks). In the message-passing paradigm, each MPI rank handles the simulation within its subdomain and communicates with the others to send data they have and to receive those necessary data that they can not compute by themselves to achieve parallelism.

4.1. Mesh Partition

For the numerical results, two different kinds of meshes have been considered: completely unstructured meshes, which have been created by using commercial software, such as ANSYS Meshing${}^{™}$ or Gambit${}^{\circledR }$, and regular unstructured meshes generated within the code using a methodology based on split Cartesian meshes. In case the output format of the external meshing software does not coincide with any of the formats admitted by the solver, we make use of the feconv mesh converter [90], to transform it. The mesh partitioning is carried out using the free METIS software package described in [59] when the unstructured grid is created with commercial software. This software provides an allotment where the number of elements of the primal mesh assigned to each MPI rank is almost the same. In this way, we achieve the best performance because computations performed by each CPU are optimized, and the number of MPI communications among CPUs is minimized. However, this method is not the most efficient when the grid is extremely fine. The most relevant reason is that the main rank has to read the entire mesh from disk and subsequently needs to perform mesh-related computations that are not related to the numerical solution but to the algorithm parallelization: generation of the staggered dual mesh topology, data to be distributed has to be prepared before being sent, the mesh has to be partitioned, and these tasks are performed by a single rank or by all ranks without parallelizing the process. Moreover, it can also yield a huge memory consumption which may not be affordable for standard nodes in HPC. To overcome this difficulty, we have implemented a parallel split Cartesian mesh generator that, as its name suggests, creates an initial Cartesian mesh of hexahedra, each of which is then divided into five tetrahedra following [91], and thus a primal mesh of tetrahedra is obtained. Figure 3 shows a cube cartesian mesh of 32 hexahedra in each direction, partitioned in 8 submeshes and a detail of the division of the hexahedra into tetrahedra. To make good use of the parallel implementation, once the total number of hexahedra and the number of partitions in each direction are set, each MPI rank meshes only its own subdomain, using the parallel grid generator previously described. For the MPI parallelization of the explicit finite volume schemes on the dual mesh and for the implicit pressure system, we need to consider a layer of shadow elements around each subdomain in order to exchange the required cell averages (FV) and vertex data (FE) between neighbor processors.

4.2. Neighbour Elements

When the mesh is partitioned, the elements $\left\{T_k\right\}$, and vertices $\left\{V_n\right\}$ of the primal mesh are distributed on the different MPI ranks, and after that, each MPI rank builds the corresponding dual mesh, and the resulting dual cells $\left\{C_i\right\}$ are subsequently distributed. To compute, for example, the velocity or density fields at cell $C_i$ or the pressure field at each primal element $T_k$, some information is required from neighbour finite volumes or finite elements, and there is no guaranty that both the element where the field is computed and the neighbour element belong to the same MPI rank. Then, it is necessary to communicate this information from one rank to the other. For this purpose, each rank will have a set of shadow primal elements and shadow dual cells, which contain data coming from neighbour processors. Each rank computes each field only in its own elements and then exchanges information with the neighbour processors that have these elements as a shadow. The exchange of information is performed with non-blocking communications based on standard calls to the MPI library. In the code, there are maps that go from local elements or vertices of the primal mesh (own and shadow elements of the partitioned mesh for each rank) to global elements in the original mesh without the partition and vice versa. In this way, it is easy to go from the original mesh elements to each partition and from the partition elements to the initial mesh. This neighbor search is completed at the beginning and is invariant for the entire simulation.

4.3. Periodic Boundaries

When periodic boundaries are considered, the distribution of the primal mesh elements $\left\{T_k\right\}$ is the same as without periodic boundaries, but the distribution of the vertices $\left\{V_n\right\}$ is slightly modified. First of all, notice that if a vertex $V_n$ is on one of the periodic boundaries, another vertex, $V_n^{\prime }$, is located in the same position on the matching periodic boundary. Both vertices are identified as the same one (merged), and the value of any field at both points must be equal. Then, if a vertex $V_n$ is located in one of the periodic boundaries and belongs to one rank, it is considered that $V_n^{\prime }$ is also own, and if the vertex $V_n$ is a shadow one, then $V_n^{\prime }$ is also shadow. In this way, although each finite element $T_k$ belongs only to one rank, the boundary vertex on a periodic boundary could belong to more than one rank.

4.4. Scalability

To illustrate the advantages of using parallelization in the family of hybrid methods, a scalability study is performed by using the Taylor-Green vortex (TGV) benchmark in a three-dimensional domain $\mathrm{\Omega }={[0,2\pi ]}^3$ (see [46,86] and Section 5.1.2 for further details about this test), with a grid created by using the split Cartesian method. First of all, a strong scaling study is considered. If the problem size n is fixed, the most common measure of scalability is the simulation speedup which can be measured as a function of the number of processors used in the simulation p. Thus, the speedup is given by

$$S\left(p\right)=\frac{T_{\mathrm{serial}}}{T_{\mathrm{parallel}}\left(p\right)},$$

(46)

where $T_{\mathrm{parallel}}\left(p\right)$ is the execution time required by p processors, and $T_{\mathrm{serial}}$ is the time necessary to execute the serial code. The linear speedup is the ideal situation where the speedup is equal to the number of processors ($S\left(p\right)=p$), and every processor contributes $100\%$ of its computational power. Moreover, the efficiency, $E\left(p\right)$, can be defined as the ratio of the speedup to the ideal linear speedup, that is,

$$E\left(p\right)=\frac{S\left(p\right)}{p}=\frac{T_{\mathrm{serial}}}{p{T}_{\mathrm{parallel}}\left(p\right)}.$$

(47)

To perform the scalability study, first of all, the left plot of Figure 4 shows the wall-clock time for solving the TGV benchmark with ${256}^3$ hexahedra, that is, with 83,886,080 tetrahedra. In this case, the wall-clock time is given as a function of the number of processors. The right plot shows the efficiency in terms of the number of processors going from 16 up to 4096. All the simulations in this section have been run on the high-performance supercomputer SuperMUC-NG of the Leibniz Rechenzentrum (LRZ) in Garching, Germany.

Figure 5 shows the strong scaling study for the parallel computation of the TGV benchmark with three different meshes: ${64}^3\times 5$ tetrahedra, in blue, ${128}^3\times 5$ tetrahedra, in green, and ${256}^3\times 5$ tetrahedra, in magenta. In the left plot, the speedup computed by using (46) is plotted in terms of the number of processors and compared with the ideal wall-clock time from 16 up to 4096 CPUs. The solid lines represent the computed speedup, and the dashed line the ideal linear speedup. The right plot shows the efficiency calculated by using (47). Without doing any code optimization, the achieved efficiency is around $60\%$ on 4096 CPU cores. As it can be observed, the finer the mesh and the greater the number of elements, the better the efficiency and speedup. By increasing the number of elements, the efficiency of the MPI communication improves, as expected.

In strong scaling studies, the problem size is fixed, but in some cases, it is interesting to study how to increase parallelism to solve problems with a bigger size. In these cases, a weak scaling study can be performed and it is used the scaled speedup given by

$$S\left(p\right)=\frac{T_{\mathrm{serial}}\left(n\left(p\right)\right)}{T_{\mathrm{parallel}}(n\left(p\right),p)},$$

(48)

where $n\left(p\right)$ is the problem size chosen so that the work per processor remains constant. In Figure 6, the scaled speedup, computed by using (48), of the TGV benchmark is shown. The problem size varies from 163,840 up to 10,485,760 elements, and the number of processors is chosen between 1 and 64 to maintain the number of elements per processor constant and equal to 163,840.

It is important to highlight that the scalability study has been performed without doing any previous code optimization. Part of our future work will consist in optimizing the code to further improve its parallel efficiency. As it can be observed in this section, when we increase the number of processors, we do not reach the ideal efficiency, and it could happen for several reasons. Increasing the number of processes sometimes will not improve the parallel code scaling but may even decrease its performance, and if this happens by increasing the number of CPUs, there is a communication overhead. The parallel algorithm should not generate communication overhead, i.e., each CPU should spend more time performing computations than sending information to achieve an efficient algorithm. And, in addition, sometimes there is a waiting overhead because ranks must wait for data from other ranks to perform a computation.

To guess where to focus our efforts on improving efficiency, we have solved the TGV benchmark with 83,886,080 elements considering the hybrid FV/FE method described in this paper, but we have also solved it considering only finite elements and only finite volumes. The left plot in Figure 7 compares the execution time when we solve with the hybrid FV/FE method with solid blue line, only with finite volumes with dashed red line and only with finite elements with dash-dotted cyan line. As expected, the wall-clock time decreases when using only FV. In the right plot of Figure 7, the efficiency comparison is shown, and, as we suspected, the efficiency when solving with the FE method alone worsens drastically (almost half that with the hybrid methodology). This study suggests that we should optimize over all the FE code. In the projection stage, the pressure equation is solved using a finite element method, and the resulting linear systems are symmetric and positive definite. For this reason, we have used a matrix-free conjugate gradient (CG) method where the parallelization of the matrix-vector product is efficient and almost straightforward. The bottleneck appears when each MPI rank should communicate in each CG iteration to the master rank the scalar norm of the residual produced, and the master checks if the tolerance imposed for the iterative method to finish is reached. Then we have a waiting overhead and an all-to-all communication of the norm of the residual that is very time-consuming on large CPU numbers, and that should be improved in the future. Despite that, we consider that the efficiency is good enough, taking into account that the mesh is really fine and that no optimization has been performed on the code so far. Future work will concern the extension of the parallel implementation of the code to a mixed MPI + OpenMP parallelization, in order to reduce the number of ranks involved in global all-to-all communications.

5. Numerical Results

This section aims at illustrating the behavior of the proposed methodology. To this end, several classical benchmarks are reproduced. The different tests carried out are organized according to the mathematical models solved, i.e., incompressible Navier-Stokes equations, weakly compressible Navier-Stokes equations, fully compressible Navier-Stokes equations, and shallow water equations.

Let us note that, for stability purposes and unless stated otherwise, the time step is dynamically computed according to a CFL-type restriction

$$\Delta t=\underset{C_i}{min}\left\{\Delta t{}_i\right\},\Delta t{}_i=\mathrm{CFL}\frac{r_i^2}{({\left|\zeta \right|}_{max}+c_{\alpha })r_i+2{\left|\nu \right|}_{max}}$$

(49)

for the incompressible and weakly compressible Navier-Stokes equations,

$$\Delta t=\underset{C_i}{min}\left\{\Delta t{}_i\right\},\Delta t{}_i=\mathrm{CFL}\frac{r_i^2}{({\left|\zeta \right|}_{max}+c_{\alpha })r_i+max\left\{2{\left|\nu \right|}_{max},\frac{c_v\gamma }{\kappa \rho }\right\}}$$

(50)

for the fully compressible Navier-Stokes equations and

$$\Delta t=\underset{C_i}{min}\left\{\Delta t{}_i\right\},\Delta t{}_i=\mathrm{CFL}\frac{r_i^2}{({\left|\zeta \right|}_{max}+c_{\alpha })r_i+{\displaystyle \frac{8}{3}}{\left|\nu \right|}_{max}}$$

(51)

for the shallow water equations. Above, we have denoted by ${\left|\zeta \right|}_{max}$ the maximum absolute eigenvalue related to the convective terms, while ${\left|\nu \right|}_{max}$ corresponds to the one of diffusion terms and is activated only if they are computed using an explicit approach. Furthermore, $c_{\alpha }$ is the artificial viscosity parameter, set to zero unless a different value is explicitly given, $c_v$ is the heat capacity at constant volume, and $r_i$ denotes the incircle diameter of each dual control volume.

The international system of units (SI) is used throughout all test cases.

5.1. Incompressible Navier-Stokes Equations

The first mathematical model considered here are the incompressible Navier-Stokes equations. First, the accuracy of the proposed hybrid FV/FE methodology is assessed using the Taylor-Green vortex benchmark in 2D. Then, we will show the ability of the scheme to solve also complex three-dimensional flows. To this end, the first test to be addressed is the three-dimensional Taylor-Green vortex which presents a complex laminar to turbulent transition behavior of the fluid flow. Several Reynolds numbers are analyzed, comparing the results obtained with the hybrid FV/FE methodology against available DNS reference data. The last test presented here is the classical lid-driven cavity benchmark in 3D, for whose verification numerical reference data are employed.

5.1.1. Taylor-Green Vortex in a Two-Dimensional Domain

To assess the accuracy of the hybrid FV/FE methodology, we consider the Taylor-Green vortex benchmark in 2D, whose known exact solution is given by

$$\mathbf{u}\left(\mathbf{x},t\right)=\left(\begin{array}{c}\hfill sin\left(x\right)cos\left(y\right)\\ \hfill -cos\left(x\right)sin\left(y\right)\end{array}\right),p\left(\mathbf{x},t\right)=\frac{1}{4}\left(cos\left(2x\right)+cos\left(2y\right)\right).$$

(52)

The computational domain $\mathrm{\Omega }={\left[0,2\pi \right]}^2$ is discretized with the mapped triangular primal meshes described in

Table 1. Description of the meshes and time steps used to perform the convergence analysis with the Taylor-Green vortex benchmark in 2D.

Mesh	Elements	Vertices	Dual Elements	$\mathbf{\Delta }\mathit{t}$
$M_1$	128	81	208	$2.5\times {10}^{-2}$
$M_2$	512	289	800	$1.25\times {10}^{-2}$
$M_3$	2048	1089	3136	$6.125\times {10}^{-3}$
$M_4$	8192	4225	12,416	$3.125\times {10}^{-3}$
$M_5$	32,768	16,641	49,408	$1.5625\times {10}^{-3}$
$M_6$	131,072	66,049	197,120	$7.8125\times {10}^{-4}$
$M_7$	524,288	263,169	787,456	$3.90625\times {10}^{-4}$
$M_8$	2,097,152	1,050,625	3,147,776	$1.953125\times {10}^{-4}$
$M_9$	8,388,608	4,198,401	12,587,008	$9.765625\times {10}^{-5}$

. Note that, for this test case, the time step has been fixed a priori according to the mesh refinement instead of setting a constant CFL value. As shown in

Table 2. Spatial $L^2$ error norms and convergence rates at time $t=0.1$ obtained using the LADER scheme for the Taylor-Green vortex benchmark in 2D.

Mesh	${\mathit{L}}_{\mathbf{\Omega }}^2\left(\mathit{p}\right)$	$\mathcal{O}\left(\mathit{p}\right)$	${\mathit{L}}_{\mathbf{\Omega }}^2\left({\mathbf{w}}_{\mathbf{v}}\right)$	$\mathcal{O}\left({\mathbf{w}}_{\mathbf{v}}\right)$
M1	4.13E−01		1.27E−01
M2	1.16E−01	1.83	3.08E−02	2.04
M3	2.99E−02	1.96	7.63E−03	2.01
M4	7.54E−03	1.98	1.90E−03	2.00
M5	1.89E−03	1.99	4.75E−04	2.00
M6	4.74E−04	2.00	1.19E−04	2.00
M7	1.19E−04	2.00	2.97E−05	2.00
M8	2.97E−05	2.00	7.42E−06	2.00
M9	7.42E−06	2.00	1.85E−06	2.00

, the sought second order accuracy is perfectly achieved using the second order local ADER scheme in the transport-diffusion stage. The simulations have been run on 16 cores of an Intel${}^{\circledR }$ Core${}^{™}$ i9-10980XE workstation.

5.1.2. Taylor-Green Vortex in a Three-Dimensional Domain

Now, the Taylor-Green vortex in a three-dimensional domain $\mathrm{\Omega }={\left[-\pi ,\pi \right]}^3$ is considered. The initial condition [47,92], is given by

$$\mathbf{u}\left(\mathbf{x},t\right)=\left(\begin{array}{c}sin\left(x\right)cos\left(y\right)cos\left(z\right)\\ -cos\left(x\right)sin\left(y\right)cos\left(z\right)\\ 0\end{array}\right),p\left(\mathbf{x},t\right)=\frac{1}{16}\left(cos\left(2x\right)+cos\left(2y\right)\right)(cos\left(2z\right)+2).$$

(53)

Variation of the viscosity parameter according to predefined Reynolds numbers, $Re\in \left\{100,200,400,800,1600\right\}$, allows the definition of a set of test cases. Despite the initial condition to be smooth, this benchmark rapidly evolves into a turbulent flow. To properly capture the complex structures generated for high Reynolds numbers, very fine grids are required. Generation of the related auxiliary mesh data in a single core might be too RAM demanding for standard HPC nodes. To overcome this issue, the parallel mesh generator described in Section 4.1 is employed. Accordingly, we provide the code with the number of divisions along each axis direction, and each CPU will generate and know only the hexahedral submesh corresponding to its computational subdomain, plus a shadow layer needed for MPI communication. Then, the associated primal mesh made of tetrahedra is built by splitting each hexahedron into five tetrahedra, as well as the related dual mesh.

Table 3. Number of spatial divisions $N_x=N_y=N_z$ and number of MPI subdomains (CPU cores) $C_x=C_y=C_z$ along each coordinate direction of the grids used for the solution of the three-dimensional Taylor-Green vortex with different Reynolds numbers.

Re Number	${\mathit{N}}_{\mathit{x}}={\mathit{N}}_{\mathit{y}}={\mathit{N}}_{\mathit{z}}$	${\mathit{C}}_{\mathit{x}}={\mathit{C}}_{\mathit{y}}={\mathit{C}}_{\mathit{z}}$
100	128	8
200	128	8
400	128	8
800	256	16
1600	512	32

details the grid resolution and the number of MPI subdomains employed in each spatial direction for each Reynolds number (Re). The total number of elements of each simulation is, therefore, $N_{\mathrm{tot}}=5·N_x·N_y·N_z$ and the total number of CPU cores and thus MPI ranks is $C_{\mathrm{tot}}=C_x·C_y·C_z$. Periodic boundary conditions are considered everywhere.

The simulations were carried out up to a final time of $t=10$ with $CFL=0.5$ and using the local ADER methodology without any limiters or artificial viscosity. The left subplots of Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 show thirteen equidistant isosurface contours for the pressure field, $p\in \left[-0.1,0.4\right]$, at the final time. The total kinetic energy dissipation rates for the considered Reynolds numbers are presented in the corresponding right subplots and compared with available DNS reference data given in [93]. An excellent agreement is observed even for the high Reynolds number test cases.

5.1.3. Lid-Driven Cavity in 3D

The lid-driven cavity test is a classical benchmark for incompressible flow solvers [48,94]. The computational domain, $\mathrm{\Omega }={[0,1]}^3$, contains an initial fluid at rest, which starts moving thanks to the displacement of the upper boundary of the domain, whose velocity is of the form $\mathbf{v}=(1,0,0)$. Homogeneous Dirichlet boundary conditions are imposed for the velocity on the bottom and lateral boundaries. Two different Reynolds numbers, $Re=400$ and $Re=1000$ are studied using the incompressible hybrid FV/FE solver with artificial viscosity coefficient $c_{\alpha }=2$ for the convection terms. The solution obtained at time $t=40$ is reported in Figure 13. The streamlines clearly show the three-dimensional effects arising in this test case. Contour plots of the velocity components at the mid-planes of the domain allow for a quantitative comparison of the obtained solution against available reference data [48,94]. Moreover, Figure 14 depicts the comparison between the approximated and the reference solution for the horizontal and vertical velocities along the vertical and horizontal 1D cuts in the middle of the domain. A good agreement is observed for both Reynolds numbers studied in this paper.

5.2. Weakly Compressible Navier-Stokes Equations

We now consider the hybrid FV/FE method applied to the non-conservative weakly compressible Navier-Stokes Equations (7)–(9). The first two test cases are focused on natural convection applications and analyze the evolution of the flow field generated by a heated bubble embedded in an initial fluid at rest. The third test case analyses the sound waves generated by the unsteady viscous laminar flow around a circular cylinder and corresponds to a medium Mach number problem.

5.2.1. Rising Bubble in 2D

As first test case addressing thermal convection-driven problems in compressible flows, we consider a two-dimensional rising bubble following classical benchmarks introduced in [30,50,95,96,97,98]. In particular, we define the computational domain $\mathrm{\Omega }=\left[0,2\right]\times \left[0,3\right]$, an initial fluid at rest and a thermal bubble centered in ${\mathbf{x}}_b=\left(1,1.5\right)$:

$$\rho \left(x,y,0\right)=\left\{\begin{array}{cc}1-{\displaystyle \frac{1}{2}}{e}^{-\frac{1}{2}{\left(\frac{r}{0.2}\right)}^2}\hfill & \hfill \mathrm{if}{r}^2\le 0.1,\\ 1-{\displaystyle \frac{1}{2}}{e}^{-\frac{5}{4}}\hfill & \hfill \mathrm{otherwise};\end{array}\right.p\left(x,y,0\right)={10}^5+\mathbf{g}·\mathbf{x},$$

(54)

where $r=\sqrt{{\left(x-1\right)}^2+{\left(y-1.5\right)}^2}$ is the radius with respect to ${\mathbf{x}}_b$. We furthermore set the viscosity $\mu ={10}^{-2}$, the thermal conductivity $\lambda =2.38\times {10}^{-2}$, the heat capacity ratio $\gamma =1.4$, the heat capacity at constant volume $c_v=717.5$, and the gravity $\mathbf{g}={\left(0,-9.81\right)}^T$. The final simulation time is $t=1.5$, and the domain has been discretized using 501,680 primal elements run on 72 CPU cores of HPC2-UNITN supercomputer. Periodic boundary conditions are defined in the y-direction, while the fixed values imposed in the upper and bottom boundaries for the velocity and temperature fields are assumed to match the initial conditions.

According to the magnitude of the physical viscosity and to reduce the computational cost of the algorithm, the viscous terms have been computed implicitly. To this end, the viscous subsystems are defined following [48,88], and continuous finite elements are used for its solution, see [88] for details. The explicit transport stage has been solved using the local ADER method without limiters and $c_{\alpha }=0.5$.

The contour plot of the temperature field obtained at times $t\in \left\{0.5,1,1.5\right\}$ is shown in Figure 15. As expected, the bubble starts rising due to gravity and natural convection effects, generating a non-zero velocity field, while the maximum temperature decreases due to heat conduction (see [50]). Let us note that, even if this test case is included in the weakly compressible flows section of this paper, the simulation has been carried out also using the hybrid FV/FE solver for fully compressible flows. A good agreement can be observed between both methods. Due to the low Mach number regime and taking into account the computational cost of both solvers, it is advisable to employ the weakly compressible method within this kind of regime. Consequently, for the three-dimensional natural convection test presented in the next subsection, we will only consider the weakly compressible solver. Further comparisons, in the context of natural convection problems, of the proposed methodology with a very high order DG numerical scheme for all Mach number flows [48], can be found in [30].

5.2.2. Rising Bubble in 3D

This test corresponds to the 3D version of the former rising bubble. We consider the computational domain $\mathrm{\Omega }=\left[0,2\right]\times \left[0,2\right]\times \left[0,5\right]$ discretized with a primal mesh made of $10,240,000$ tetrahedral elements built using the parallel split Cartesian grid generator with $N_x=N_y=160$, $N_z=400$. The initial conditions correspond to (54) with $r^2={(x-1)}^2+{(y-1)}^2+{(z-1)}^2$. To complete the setup of this test, we define $\mu ={10}^{-3}$, $\lambda =2.38\times {10}^{-2}$, $\gamma =1.4$, $c_v=717.5$, $\mathbf{g}={\left(0,0,-9.81\right)}^T$. Periodic boundary conditions have been set in x and y directions, while Dirichlet boundary conditions for the conserved variables and Neumann boundary conditions for the pressure are considered in the upper and bottom boundaries. The simulation has been run on the high-performance supercomputer SuperMUC-NG using 2000 CPU cores. Figure 16 depicts the thermal bubble rising up in the domain for times $t\in \left\{0,1,2,3,4,5\right\}$. Due to the differences of temperature between the bubble and the remaining computational domain, a velocity field conveying the heated flow upwards is generated, see Figure 17. Accordingly the initially round bubble rises and gradually deforms into a mushroom-shaped structure. Despite no reference solution is available for this test case, a qualitative comparison against numerical and experimental data can be performed, and the symmetry of the bubble can be checked, confirming the correct behavior of the applied methodology, see (e.g., [50]) for comparison.

5.2.3. Aeolian Tones Generated by the Compressible Flow around a Circular Cylinder

It is well known that the unsteady flow of a von Karman vortex street in the wake of a circular cylinder generates aeolian tones within a compressible fluid, see (e.g., [99]) and references therein. This section reproduces the test according to the setup detailed in [99,100]. Our computational domain $\mathrm{\Omega }$ is a circle of radius $R=300$ centered in the origin, and the embedded circular cylinder with the same center has a diameter of $d=1$. The Reynolds number based on the diameter of the cylinder is $Re=\frac{{\rho }_{\infty }{\mathbf{v}}_{\infty }d}{\mu }=150$, while the Mach number of the flow is set to $M_{\infty }=0.2$, hence it is sufficient to solve the weakly compressible Navier-Stokes equations. Heat conduction is neglected in our test. The ratio of specific heats is set to the usual value of $\gamma =1.4$, the freestream density, pressure, and velocity are ${\rho }_{\infty }=1$, $p_{\infty }=1/\gamma$, and $v_{\infty }=0.2$, respectively.

The simulation is run on an unstructured triangular mesh composed of 667,340 triangles with a characteristic mesh spacing of $h=0.05$ at the cylinder and $h=4$ at the outer boundary. As already mentioned before, the computational domain is a circle with radius $R=300$, and the outer region $200\le \parallel \mathbf{x}\parallel \le 300$ is used as a sponge layer in order to allow the sound waves to leave the domain without spurious reflections. Within the sponge layer, the degrees of freedom of the pressure ${\widehat{\mathrm{\Pi }}}_i$ are simply overwritten in each timestep as follows:

$${{\tilde{\mathrm{\Pi }}}_i}^{n+1}=(1-m){{\widehat{\mathrm{\Pi }}}_i}^{n+1}+m{p}_{\infty },$$

(55)

where ${{\widehat{\mathrm{\Pi }}}_i}^{n+1}$ are the degrees of freedom (DoF) obtained after the projection stage of the hybrid FV/FE scheme, and ${{\tilde{\mathrm{\Pi }}}_i}^{n+1}$ is the final value assigned to each DoF after the application of the sponge layer. The absorption coefficient is set to $m=3\times {10}^{-5}$. The simulation is performed on 2400 CPU cores of the SuperMUC-NG supercomputer until $t=900$.

The sound pressure field generated by the von Karman vortex street is depicted in Figure 18 and agrees qualitatively very well with the sound field obtained in [99,100]. The time series of the velocity component $v_2({\mathbf{x}}_p,t)$ in the point ${\mathbf{x}}_p=(10,0)$ is shown in Figure 19. From this time signal, we obtain a Strouhal number of $St=\frac{fd}{u_{\infty }}=0.179$, which can be compared with the reference value of $St=0.183$ found by Müller in [99] and the other reference values given therein. The Strouhal number obtained in [100] was $St=0.182$. The error with respect to [99,100] is of the order of $2\%$, hence we obtain a rather good quantitative agreement with reference results published in the literature.

5.3. Fully Compressible Navier-Stokes Equations

This section shows some computational results obtained with the hybrid FV/FE method applied to the fully compressible Navier-Stokes equations in two and three space dimensions, see also [31], where this method was introduced for the first time. The first test is a two-dimensional Kelvin-Helmholtz instability run at a low Mach number, while the second test consists of two shock tube problems in order to show that the proposed scheme is indeed an all Mach number flow solver that is able to solve both low Mach number flows as well as high Mach number flows with shocks and discontinuities.

5.3.1. Kelvin-Helmholtz Instability

The computational domain used for this simulation is the square $\mathrm{\Omega }=[0,1]\times [0,1]$, with periodic boundary conditions in the $x-$direction and slip wall boundaries in the $y-$direction. The initial condition of this test problem is given by

$$\rho =1-\frac{1}{4}tanh\left(25y\right),v_1=-\frac{1}{2}tanh\left(25y\right),v_2=\frac{1}{100}sin\left(2\pi x\right)sin\left(\pi y\right),p=\frac{{10}^4}{\gamma }.$$

(56)

The average Mach number of this test problem is, therefore, $M=0.01$. The viscosity coefficient is set to $\mu =0.0001$, while thermal conduction is neglected, i.e., $\lambda =0$. A second order accurate local ADER scheme without limiter is employed to discretize the nonlinear convective terms. The simulation is run on 2400 CPU cores of the SuperMUC-NG supercomputer until a final time of $t=8$, using an unstructured primal simplex mesh with characteristic mesh spacing $h=1/1000$ composed of a total of 2,279,386 triangular elements. The time evolution of the density contours is shown in Figure 20. One can notice the development of the typical vortex structures that are expected for this type of physically unstable shear flow.

5.3.2. Riemann Problems

In this section, we solve two Riemann problems that are well-known standard benchmark problems for numerical methods applied to nonlinear hyperbolic conservation laws, but here we are using an unstructured tetrahedral mesh in three space dimensions in order to check our hybrid FV/FE algorithm in 3D. The computational domain is $\mathrm{\Omega }=[-0.5,0.5]\times {[0.0125,0.0125]}^2$ with Dirichlet boundary conditions in the $x$-direction and periodic boundaries in the $y$- and $z$-directions. The unstructured mesh is a regular one, with $5\times 512\times 32\times 32$ elements. The first test problem is the classical Sod shock tube [101], while the second test problem is taken from the textbook of Toro [23] and corresponds to Riemann problem number 4. The initial condition for a generic quantity q takes the form

$$q(x,0)=\left\{\begin{array}{cc}{q}_L\hfill & \mathrm{if}x\le x_c,\\ q_R\hfill & \mathrm{if}x>x_c.\end{array}\right.$$

(57)

According to the above-mentioned references, the left and right states for density, velocity component $v_1=u$ and pressure are briefly summarized in

Table 4. Initial states left and right, initial position of the discontinuity $x_c$, and final simulation time t for each Riemann problem.

Test	${\mathit{\rho }}_{\mathit{L}}$	${\mathit{\rho }}_{\mathit{R}}$	${\mathit{u}}_{\mathit{L}}$	${\mathit{u}}_{\mathit{R}}$	${\mathit{p}}_{\mathit{L}}$	${\mathit{p}}_{\mathit{R}}$	${\mathit{x}}_{\mathit{c}}$	t
RP1	1	$0.125$	0	0	1	$0.1$	0	$0.2$
RP4	$5.99924$	$5.99242$	$19.5975$	$-6.19633$	$460.894$	$46.095$	$-0.2$	$0.035$

, including the final time of the simulation $t_{\mathrm{end}}$ and the initial position of the discontinuity $x_c$. For the 3D setup we set the tangential velocity components to zero, i.e., $v_2=v_3=0$. For this test we use an inviscid fluid, setting $\mu =\lambda =0$ and $\gamma =1.4$. Simulations are carried out on 16 CPU cores of an Intel${}^{\circledR }$ Core${}^{™}$ i9-10980XE workstation with 64 GB of RAM.

In Figure 21, the numerical results obtained with our 3D hybrid FV/FE method are compared with the exact solution of the Riemann problem provided in [23]. We obtain a very good agreement between the exact solution and numerical results for all test problems and all flow variables.

5.4. Shallow Water Equations

The last model discretized using the hybrid FV/FE methodology corresponds to the Saint-Venant equations. We consider two different test cases of dambreaks with known numerical and experimental reference solutions, thus allowing for a careful validation of the proposed algorithm.

5.4.1. 3D Dambreak on a Dry Plane

We consider the three-dimensional dam break on a dry plane proposed in [102], where experimental data and numerical results are provided. The physical domain consists of a tank of dimensions $1\times 2\times 0.7$ connected to a flat plate through a removable gate that can be opened in less than $0.1$ s, simulating an instantaneous dam break, see Figure 22 for further details on the computational domain. The rupture of the gate leads to a water front traveling over the dry plane, while in the meantime, a rarefaction wave propagates backward inside the tank. In the numerical framework, the test can be reproduced by considering symmetry boundary conditions for the tank walls and free outflow in the three remaining boundaries delimiting the initially dry plane. Within the model, we set the friction coefficient to $n={10}^{-3}$ so that experimental data can be correctly reproduced. Regarding the algorithm, the generation of spurious oscillations is avoided by using a Rusanov-type dissipation in the FE discretization, see [74] for further details.

The simulation has been run up to time $t=2$ in a primal mesh made of 2,178,836 simplex elements using 2400 CPU cores of SuperMUC-NG. Figure 23 depicts the surface elevation obtained at times $t\in \left\{0.1,0.2,0.3,0.5,1.0,2.0\right\}$. Moreover, the time evolution for the wave gauges described in

Table 5. Coordinates of the wave gauges of the three-dimensional dambreak on a dry plane for a bi-dimensional domain $\mathrm{\Omega }=[-1,1]\times [-1,2]$.

Wave Gauge	−5 A	−3 A	−2 A	0	1 A	8 A
x	$-0.82$	$-0.62$	$-0.42$	0	$0.02$	$0.722$
y	0	0	0	0	0	0

is shown in Figure 24 jointly with the experimental data provided in [102]. To ease the comparison with numerical data from different numerical approaches, we also include the solution obtained in [103] using a three-dimensional smooth particle hydrodynamics (SPH) scheme for the weakly compressible free surface Navier-Stokes equations, the numerical results for the diffuse interface method (DIM) for the weakly compressible Euler equations proposed in [104] and the approximation obtained in [102] with a Godunov-type finite volume scheme for the shallow water equations. As expected, the hydrostatic Saint Venant model provides worse results than the non-hydrostatic fully 3D Navier-Stokes solvers at short times due to the hydrostatic pressure assumption in the shallow water equations, but for larger times, the results greatly improve while presenting a reduced computational cost.

Note that this test case has already been run with the hybrid FV/FE scheme using a coarse mesh in [74], assessing the capability of the hybrid FV/FE methodology to address practical applications. Nevertheless, some small structures could not be captured with the grid resolution employed in [74]. To this end, a much finer grid and thus a larger number of parallel processors is needed. The new results presented in this paper confirm the ability of the hybrid methodology to accurately solve dambreak problems showing an excellent agreement with experimental and other numerical reference data, thus demonstrating the importance of the massive parallel computing to attain the level of accuracy sought in practical applications.

5.4.2. Dambreak into a Channel with a 45${}^{\circ }$ Bend (CADAM)

The second dambreak problem studied has been first put forward in [105] in the framework of the European Concerted Action on Dam Break Modeling (CADAM). Figure 25 sketches the geometry considered in the physical experiment, consisting of a reservoir connected to a rectangular channel with a ${45}^{\circ }$ bend. We assume the origin of the coordinate system to be placed at the lower left corner of the reservoir, with dimensions ${\mathrm{\Omega }}_{reservoir}=[0,2.39]\times [0,2.44]\times [-0.33,0.4]$, while the bottom of the first channel section is located at $z=0$ so that a step of height $0.33$ meters is placed at the gate between both subdomains. Like for the previous test case, a physical gate is initially placed at the reservoir exit separating it from the dry channel. Initially, the tank is filled with water at rest up to a height of $h=0.58$ m from the bottom, that is, with the free surface being $\eta =0.25$ m. The gate is then removed, simulating the rupture of a dam. Consequently, the water flows along the channel, hitting the bend walls and yielding to oblique waves propagating down the channel as well as a backward front traveling towards the reservoir. Further details on the experimental test description can be found in [105].

For the setup of the numerical test case, we consider symmetry boundary conditions at the walls and an outflow boundary condition at the end of the channel. Following the guidelines in [105], the friction coefficient is set to $n={10}^{-3}$. The primal mesh used is composed of 2,495,116 triangular elements to be divided into 2400 CPU cores. The free surface elevation for $t\in \left\{0.5,1.0,2.5,4.0,6.0,8.0\right\}$ is represented in Figure 26. All the structures presented are in good qualitative agreement with the experiments in [105], the 3D SPH simulations [106], and the solution obtained for the 3D third-order one-step WENO finite volume scheme [104]. The MPI mesh partition is also depicted in the last subfigure of Figure 26. To better analyze the numerical solution of the hybrid methodology and to ease its comparison with available reference data, we also report the time evolution of the water surface at the fixed locations given in

Table 6. Coordinates of the gauges G1–G8 of the CADAM test case.

Gauge	G1	G3	G4	G5	G6	G7	G8
x$\left[m\right]$	1.59	4.24	5.74	6.74	6.65	6.56	7.07
y$\left[m\right]$	0.69	0.69	0.69	0.72	0.80	0.89	1.22

and already depicted in Figure 25, see Figure 27. A good agreement is observed for all wave gauges. The solution reasonably matches the expected discharge of the reservoir, and the arrival times for the main wave front are well captured. Moreover, the hybrid FV/FE method is able to reproduce the backward traveling wave in G5 and perfectly matches the experimental time the main wave front arrives at G4, thus performing better than the available 3D solutions of DIM and SPH methodologies and than the pure Godunov-type FV scheme for the shallow water equations.

6. Conclusions

This paper presents a massively parallel implementation of a novel and versatile family of pressure-based staggered semi-implicit hybrid finite volume/finite element schemes for computational fluid dynamics using the MPI standard. The proposed hybrid FV/FE scheme can be applied to a wide range of incompressible flows [28,29], weakly compressible low Mach number flows [30], compressible all Mach number flows [31] as well as to the shallow water equations describing geophysical free surface flows [74].

The parallelization of the scheme, as well as the distribution of data and workload, is carried out using MPI in combination with spatial domain decomposition. For general unstructured meshes, we use the standard METIS software package [59] to distribute the spatial domain onto the various MPI ranks. While MPI parallelization via domain decomposition is straightforward for the explicit finite volume part of the hybrid FV/FE scheme, special care is needed concerning the efficient parallel solution of the implicit pressure system via finite elements. Since the linear algebraic systems that are generated by the finite element discretization of the pressure equation are always symmetric and positive definite, throughout this paper, we make use of a matrix-free implementation of the conjugate gradient (CG) method, which allows an efficient and almost straightforward parallel implementation of the matrix-vector products that are needed by the matrix-free CG method. The structure of the resulting parallel implicit solver is in each CG iteration similar to one of a parallel explicit scheme, with MPI data transfer via a layer of shadow/halo zones around each MPI rank. The only bottleneck is the global communication needed for the exchange of the scalar norm of the residual in the CG algorithm, but which is the least global communication effort that one can imagine in implicit pressure-based solvers. Following this parallel programming paradigm, the code can be run on modern massively parallel distributed memory supercomputers. To validate the efficiency of the parallel implementation, a strong scaling test has been performed for the three-dimensional Taylor-Green vortex problem up to 4096 CPU cores. For sufficiently large problem sizes, the total MPI efficiency was still about 60% even on 4096 CPU cores for the complete semi-implicit hybrid FV/FE algorithm. The largest simulation shown in this paper was the 3D Taylor-Green vortex, run on 32,768 CPU cores and 671 million tetrahedral elements.

The algorithm has been tested on a wide range of test cases, ranging from incompressible flows over weakly compressible low Mach number flows up to compressible flows with shock waves and free-surface flows described by the shallow water equations.

Future work will concern the extension of the present parallel semi-implicit hybrid FV/FE method to more complex PDE systems, like the MHD equations [107], the GPR model of continuum mechanics [108,109,110], as well as to first order hyperbolic reformulations of nonlinear dispersive systems [111,112,113]. We will also consider fluid-structure interaction problems (FSI) on moving meshes via a direct Arbitrary-Lagrangian-Eulerian (ALE) approach on moving conforming simplex meshes [114,115].