# Novel adaptive finite volume method on unstructured meshes for time-domain wave scattering and diffraction

Tarek Ghoudi, M. Shadi Mohamed, Mohammed Seaid

► To cite this version:

Tarek Ghoudi, M. Shadi Mohamed, Mohammed Seaid. Novel adaptive finite volume method on unstructured meshes for time-domain wave scattering and diffraction. Computers & Mathematics with Applications, 2023, 141, pp.54-66. 10.1016/j.camwa.2023.03.025. hal-04546410

HAL Id: hal-04546410

https://hal.science/hal-04546410v1

Submitted on 17 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Novel adaptive finite volume method on unstructured meshes for time-domain wave scattering and diffraction

Tarek Ghoudi<sup>a</sup>, M. Shadi Mohamed<sup>b,*</sup>, Mohammed Seaid<sup>c</sup>

## ARTICLE INFO

### Keywords:

Wave equation

Finite volume method

Vertex-centered discretization

Wave scattering

Wave diffraction

Mesh adaptation

## ABSTRACT

A new adaptive finite volume method is proposed for the simulation of the wave problems in the time domain. The transient wave equations are discretized in time and space. A vertex-centered finite volume method is constructed with both cell-centered and edge-midpoint of each control volume. We then propose a mesh adaptation procedure based on energy-norm error-estimates, which significantly increases the efficiency of the method. The proposed approach is accurate in capturing the details of the scattered and diffracted waves with highly refined elements following the evolution of the wave patterns. This is a critical feature of the approach as waves can propagate in vast domains and we successfully refine the mesh only where needed. Unlike many other methods for evolution problems in which the differential operator is solved after each error estimation, the proposed approach allows for multiple adaptations of meshes within a single error estimation. This nested adaptive finite volume method requires only treatment of conformity in meshes for multiple adaptations which is dealt with using a Newest-Vertex-Bisection algorithm. Comparisons with the finite element and with reference solutions are considered for progressive radial waves, waves reflection and a wave scattering and diffraction around barriers. The proposed approach results are more efficient and highly accurate, hence, the significant potential when applied to the time-domain simulation of wave problems.

## 1. Introduction

Dealing with subsurface imaging and seismic modelling as well as the related inverse problems often requires recovering the wave field predicted by the wave equation on large computational domains and may involve complicated geometries. A few examples can incorporate earth curvature in long-offset crustal seismic surveys [39], the full-waveform inversion from irregular surfaces [7,37] and modelling waveforms from deviated boreholes [45]. In the vast majority of such applications, the wavefield can only be recovered using numerical techniques. Numerical solutions of the wave equation can be achieved using various numerical methods, such as the finite element method (FEM) [44,23,28], the finite volume method (FVM) [49] and the finite difference method [40,6]. Out of these methods, the FEM is well known for its flexibility in dealing with complex geometries and non-homogeneous material properties. In the present work, we argue that the FVM which to an extent is a less well-known technique for solving the wave equation, can be as flexible as the FEM for solving wave problems in the

time domain [13]. Furthermore, we show that the FVM can bring other advantages into the numerical solution of the wave equation. The FEM is known to suffer from dispersion error when solving wave problems [12]. Missing to capture the accurate wavelength does not only reflect on the fine scale details of the solution, but it also contaminates the global features of the numerical solution [12]. This behavior is intensively studied in the literature and is known by the pollution error [27]. The pollution cannot be avoided in problems with multi-spatial dimensions, however, it can be reduced by either increasing the discretization resolution or introducing special measures into the FEM [3,14]. Indeed, using isogeometric analysis within the finite element framework improves the pollution behavior, see for example [30]. This has motivated different researchers to investigate multiple aspects of the method such as the iterative solution of the linear system [16,15] and the boundary condition implementation [31] and adaptivity [43]. The improvement in the pollution has also been reported in the context of the boundary element method [10]. However, all standard computational methods are to some extent affected by the pollution error [25]. In this paper

<footnote><sup>a</sup> CMLA, ENS Paris-Saclay, CNRS, Université Paris-Saclay, 94235, Cachan, France</footnote>

<footnote><sup>b</sup> School of Energy, Geoscience, Infrastructure and Society, Heriot-Watt University, Edinburgh EH14 4AS, UK</footnote>

<footnote><sup>c</sup> Department of Engineering, University of Durham, South Road, Durham DH1 3LE, UK</footnote>

<footnote>* Corresponding author.
E-mail addresses: ghoudi@math.univ-paris13.fr (T. Ghoudi), m.s.mohamed@hw.ac.uk (M.S. Mohamed), m.seaid@durham.ac.uk (M. Seaid).</footnote>
we show that the FVM is less affected by the pollution errors compared to the FEM. In addition, to improve the FVM further we propose an adaptive approach. The proposed approach is then studied for wave applications with special attention paid to wave scattering and diffraction.

The finite volume method can be efficiently used for solving time-dependent problems [47,22]. The FVM is used to approximate the wave-field as a set of cell/element averaged values [33]. The local values within one cell are interpolated from these averaged values. High accuracy spatial interpolation can be achieved using high-order polynomials [20,19] while, for the time discretization, the Runge-Kutta schemes are often used [33]. It is shown that high accuracy can be achieved using the FVM on unstructured meshes when solving the wave problems [33,20]. High-order polynomial reconstructions are successfully implemented to deal with heterogeneous problems in two- and three-dimensional problems [49,19] and used to model vibro-acoustic applications [48]. The two-dimensional FVM has also been successfully implemented with the Roe linearization to solve nonlinear wave diffraction problems [42]. To quantify the uncertainty of wave propagation in a heterogeneous media the FVM is coupled with the multi-level Monte Carlo method [35]. Instead of transforming the wave equation to a first-order hyperbolic system, an upwind scheme is proposed based on a Riemann-type problem and it can directly deal with the second order system [4]. Although this has been initially implemented using the finite difference method but it can be equally used for finite volume approximations. Notice that the wave problems investigated in the present study should not be confused with wave problems governed by hyperbolic systems of conservation laws for which well-established Riemann-based solvers. Here, the hyperbolic nature of the problem is coming from the second-order derivative in time but the spatial derivative is still governed by the elliptic part. Furthermore, Riemann-based solvers require in their formulation, the eigen-structures associated with the problem to construct the numerical fluxes which are not available in the considered problem. The wave problem in this study has been widely solved using finite element methods and our aim is to propose an adaptive finite volume method solving the wave problem which conserves all those physical properties required when it is solved in conjunction with hyperbolic systems of conservation laws for mass, momentum and energy.

The efficiency of time integration schemes can be improved as the Butcher barrier [46]. This can be achieved by exploiting the arbitrary high-order derivatives approach [41] combined with the FVM. The resulting method has been successfully applied to wave problems on unstructured meshes [20,49,11]. This class of finite volume methods not only have good conservation characteristics but also can deal with complex grids. However, these methods have been developed for solving first-order hyperbolic systems of conservation laws and mainly implemented on structured meshes using cell-centered reconstructions [13]. The purpose of the current study is to propose a vertex-centered finite volume method for solving wave problems for which the governing equations involve second-order differential operators in space and time. Here, the wave equation is discretized in space using a vertex-centered method involving both cell-centered and edge-midpoint of each control volume in its reconstruction on unstructured meshes. For the time integration of the associated semi-discrete system we apply the implicit Newmark scheme [36]. Furthermore, the novel contribution of this work is the implementation of a multilevel mesh adaption procedure in the proposed vertex-centered finite volume method. Here, we consider the flux recovery technique and error estimates in the energy norm as those investigated in [34,8] among others. Clearly, changing the quantity of interest considered in the error estimate can result into meaningful changes in the resulting mesh refinement [21]. It has also been reported in the literature that mesh labeling can have an impact on the refinement but in the context of fracture mechanics [29]. In the view of the a posteriori error estimators, a multilevel adaptive algorithm is presented for solving transient wave problems accurately and

efficiently. Unlike many methods for evolution problems in which the differential operator is solved after each error estimation, the proposed finite volume method allows for multiple adaptations of meshes within a single error estimation. This nested adaptive finite volume method uses techniques investigated in [5,32], requires only treatment of conformity in meshes for multiple adaptations which is dealt with using the Newest-Vertex-Bisection algorithm [9]. To the knowledge of the authors, numerical assessment of adaptive finite volume methods for time-dependent wave problems is reported for the first time.

Numerical results are presented for three problems of transient wave problems including a progressive radial wave problem with known analytical solution, a damped wave reflection and a wave scattering and diffraction around embedded barriers. The considered mesh adaption criteria accurately and efficiently approximates the solutions of these problems. The obtained numerical results demonstrate good resolution of the proposed vertex-centered finite volume method with high accuracy in smooth regions and without nonphysical oscillations near the steep gradients or extensive numerical dissipation. Further, the computed results are highly accurate, hence the significant method potential when applied to the directions for the time-domain simulation of wave propagation. The problems considered here are representative examples for many practical time-dependent wave applications related to scattering and diffraction. Further, the results are compared with calculations using the finite element method on fixed meshes. The remainder of this paper is organized as follows. Formulation of the vertex-centered finite volume for transient wave problems is presented in section 2 and the mesh adaptation procedure in section 3. The implementation of the Newmark scheme for the time integration is discussed in section 4. In section 5, we examine the numerical performance of the proposed method using several test examples of wave propagation problems. The proposed method provides both high accuracy and efficiency. Final concluding remarks are summarized in section 6.

## 2. Vertex-centered finite volume method for space discretization

We present a linear wave equation in two space dimensions which would be the basis of all the analysis presented in this paper. We define an initial boundary-value problem in $\Omega \subset \mathbb{R}^2$ being an open bounded domain with Lipschitz continuous boundary $\Gamma$ evolving in $[0, T]$ which is the time interval for the wave propagation. The problem is defined as

$$ \frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - \nabla^2 u = f(t, \mathbf{x}), \quad (t, \mathbf{x}) \in [0, T[\times\Omega, \qquad (2.1a) $$

$$ \frac{\partial u}{\partial \bar{v}} + hu = g(t, \mathbf{x}), \quad (t, \mathbf{x}) \in [0, T[\times\Gamma, \qquad (2.1b) $$

$$ u(0, \mathbf{x}) = u_0(\mathbf{x}), \quad \mathbf{x} \in \Omega, \qquad (2.1c) $$

$$ \frac{\partial u}{\partial t}(0, \mathbf{x}) = v_0(\mathbf{x}), \quad \mathbf{x} \in \Omega, \qquad (2.1d) $$

where $\mathbf{x} = (\mathbf{x}, \mathbf{y})^\top$ are the Cartesian coordinates, $t$ is the time variable, $\bar{v}$ the outward unit normal on $\Gamma$ and $\mathbf{u}$ the magnitude of the transverse electric field in the direction perpendicular to the plane of numerical domain while $c$ and $h$ are given constants. The functions $f(t, \mathbf{x})$ and $g(t, \mathbf{x})$ in (2.1) are the prescribed source and the boundary functions, respectively. The functions $u_0(\mathbf{x})$ and $v_0(\mathbf{x})$ denote the given initial conditions. This model can be used to represent various linear electromagnetic and acoustic wave propagation problems [38]. For instance, applied to the scalar field in a transverse mode of electromagnetic wave propagation, it can represent an accurate and efficient solution for a short pulse propagating over long distances, see for example [17,18] and further references are therein.

For the spatial discretization of (2.1), we consider the primal and dual meshes $\mathcal{T}_h$ and $\mathcal{D}_h$ shown in Fig. 2.1, elements of $\mathcal{D}_h$ are denoted by $D_i$. For $D_i \in \mathcal{D}_h$, $d_{D_i}$ denotes its diameter. For each primal element $T_i$ we denote its three vertices by $t$, $j$ and $k$, its barycenter by $G$ and its edges $[jk]$, $[ik]$ and $[ij]$, respectively, by $\Sigma_i$, $\Sigma_j$ and $\Sigma_k$. The centers of these edges are denoted by $m_{[ik]}$, $m_{[jk]}$ and $m_{[ij]}$ and the associated outward normals are $\bar{n}_i$, $\bar{n}_j$ and $\bar{n}_k$. We denote by $V_i$ the set
<figure><img src="image_3.png" /><figcaption>**Fig. 2.1.** Illustration of the primal mesh $T_i$ (solid line) and dual mesh $D_h$ (dashed line) along with the notation used in the finite volume discretization.</figcaption></figure>

of triangles sharing the vertex $i$. We also denote by $\sigma_{ij}$, $\sigma_{ik}$ and $\sigma_{jk}$ the edges $[Gm_{[ij]}]$, $[Gm_{[ik]}]$ and $[Gm_{[jk]}]$, respectively. The associated outward normals to these edges are denoted by $\vec{n}_{ij}$, $\vec{n}_{ik}$ and $\vec{n}_{jk}$. We also use $|T_i|$ and $|\sigma_{ij}|$ to denote the area of $T_i$ and the length of $\sigma_{ij}$, respectively. Hence, integrating the first equation in (2.1) over the dual element $D_i$ we obtain

$$ \frac{1}{c^2} \int_{D_i} \frac{\partial^2 u}{\partial t^2} d\mathbf{x} - \int_{D_i} \nabla^2 u \, d\mathbf{x} = \int_{D_i} f(\mathbf{x}, t) \, d\mathbf{x}. \quad (2.2) $$

To approximate the diffusive fluxes in (2.2), we apply the divergence theorem as

$$ \int_{D_i} \nabla^2 u \, d\mathbf{x} = - \oint_{\partial D_i} \nabla u \cdot \vec{n} \, d\sigma, \\ = - \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \oint_{\partial D_j \cap T_i} \nabla u \cdot \vec{n} \, d\sigma. $$

Thus, the gradient $\nabla u$ is calculated at the interface $\partial D_i \cap T_i$ as

$$ \begin{aligned} - \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_i \in V_i}} \oint_{\partial D_j \cap T_i} \nabla u \cdot \vec{n} \, d\sigma &= - \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \nabla u \cdot \left( \oint_{\partial D_j \cap T_i} \vec{n} \, d\sigma \right), \\ &= - \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \nabla u \cdot \left( |\sigma_{ij}| \vec{n}_{ij} + |\sigma_{ik}| \vec{n}_{ik} \right), \\ &= \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \frac{1}{2} \nabla u \cdot \left( |\Sigma_k| \vec{n}_k + |\Sigma_j| \vec{n}_j \right), \\ &= -\frac{1}{2} \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \nabla u \cdot \left( |\Sigma_i| \vec{n}_i \right). \end{aligned} \quad (2.3) $$

Using the fact that the solution $u$ is assumed to be piecewise linear ($P_1$ in each control volume $T_i$),

$$ \begin{aligned} \nabla u \bigg|_{\partial D_i \cap T_i} &\approx \frac{1}{|T_i|} \left( \left( \frac{u_i + u_j}{2} \right) |\Sigma_k| \vec{n}_k + \left( \frac{u_j + u_k}{2} \right) |\Sigma_i| \vec{n}_i \right. \\ &\qquad \left. + \left( \frac{u_i + u_k}{2} \right) |\Sigma_j| \vec{n}_j \right), \\ &\approx \frac{1}{2|T_i|} \left( |\Sigma_i| \vec{n}_i u_i + |\Sigma_j| \vec{n}_j u_j + |\Sigma_k| \vec{n}_k u_k \right). \end{aligned} $$

Hence,

$$ - \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \oint_{\partial D_j \cap T_i} \nabla u \cdot \vec{n} \, d\sigma $$

<figure><img src="image_4.png" /></figure>

$$ \approx \sum_{\substack{T_j \cap D_i \neq \emptyset \\ T_j \in V_i}} \frac{1}{4|T_j|} \left( |\Sigma_i| \vec{n}_i u_i + |\Sigma_j| \vec{n}_j u_j + |\Sigma_k| \vec{n}_k u_k \right) \cdot |\Sigma_i| \vec{n}_i. $$

In each element $T_i$ with vertices $i$, $j$ and $k$, we denote by $(x_i, y_i)$, $(x_j, y_j)$ and $(x_k, y_k)$ the coordinates of the vertices $i$, $j$ and $k$, respectively. In what follows we provide details on the calculation of the gradient in (2.3) by assuming that the solution $u$ is a polynomial function of first degree in the element $T_i$ which takes the values $u_i$, $u_j$ and $u_k$ at the vertices $i$, $j$ and $k$, respectively. Thus, the solution $u$ in the element $T_i$ can be formulated as

$$ u(x, y) = \lambda x + \mu y + \xi, $$

with

$$ u_i = \lambda x_i + \mu y_i + \xi, \quad u_j = \lambda x_j + \mu y_j + \xi, \quad u_k = \lambda x_k + \mu y_k + \xi. $$

Note that this is a system of three equations with three unknowns $\lambda$, $\mu$ and $\xi$. This system can be solved analytically and the gradient of $u$ can also be calculated as

$$ \frac{\partial u}{\partial x} = \lambda = -\frac{u_i}{\psi}(y_j - y_k) + \frac{u_j}{\psi}(y_i - y_k) - \frac{u_k}{\psi}(y_i - y_j), $$

$$ \frac{\partial u}{\partial y} = \mu = u_i \left( \frac{\psi + (x_i - x_j)(y_j - y_k)}{\psi(y_i - y_j)} \right) - u_j \left( \frac{\psi + (y_i - y_k)(x_i - x_j)}{\psi(y_i - y_j)} \right) + u_k \left( \frac{x_i - x_j}{\psi} \right), $$

where

$$ \psi = (x_j - x_k)(y_i - y_j) - (x_i - x_j)(y_i - y_k). $$

Note that for the case with $y_i = y_j$ we use

$$ \mu = u_i \frac{(x_i - x_j)}{\psi} + u_j \left( \frac{\psi - (y_i - y_k)(x_j - x_k)}{\psi(y_j - y_k)} \right) + u_k \left( \frac{-\psi + (y_i - y_j)(x_j - x_k)}{\psi(y_j - y_k)} \right). $$

Hence, the gradient $\nabla u$ is approximated as

$$ \nabla u \bigg|_{\partial D_i \cap T_i} \approx \left( \phi_i u_i + \phi_j u_j + \phi_k u_k \right), \quad (2.4) $$

with

$$ \phi_i = - \frac{y_j - y_k}{\psi}, \quad \phi_j = \frac{y_i - y_k}{\psi}, \quad \phi_k = - \frac{(y_i - y_j)}{\psi}, $$

and
$$
\begin{align*}
\chi_i &= \begin{cases} \frac{\psi + (x_i - x_j)(y_j - y_k)}{\psi(y_i - y_j)}, & \text{if } y_i \neq y_j, \\ \frac{(x_i - x_j)}{\psi}, & \text{if } y_i = y_j, \end{cases} \\
\chi_j &= \begin{cases} \frac{\psi + (x_i - x_j)(y_i - y_k)}{\psi(y_i - y_j)}, & \text{if } y_i \neq y_j, \\ \frac{\psi - (x_j - x_k)(y_i - y_k)}{\psi(y_j - y_k)}, & \text{if } y_i = y_j, \end{cases} \\
\chi_k &= \begin{cases} \frac{(x_i - x_j)}{\psi}, & \text{if } y_i \neq y_j, \\ \frac{-\psi + (x_j - x_k)(y_i - y_j)}{\psi(y_j - y_k)}, & \text{if } y_i = y_j. \end{cases}
\end{align*}
$$

Note that since the solution $u$ is $P_1$ in the element $T_i$, its gradient $\nabla u$
is therefore constant in $T_i$. Thus, in each element $T_i$ the equation (2.3)
can be reformulated as

$$
\sum_{T_i \cap D_i \neq \emptyset} \oint_{V_i} \nabla u \cdot \vec{n}_{ij} \, d\sigma + \sum_{T_i \cap D_i \neq \emptyset} \oint_{V_i} \nabla u \cdot \vec{n}_{ik} \, d\sigma \\
= \sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} \left(\begin{array}{c} \lambda \\ \mu \end{array}\right) \cdot \left( |Gm_{[ij]}| \vec{n}_{ij} + |Gm_{[ik]}| \vec{n}_{ik} \right).
$$

Using the gradient approximation (2.4), the above equation becomes

$$
- \sum_{\substack{T_i \cap D_i \neq \emptyset \\ T_i \in V_i}} \oint \nabla u \cdot \vec{n}_{ij} \, d\sigma - \sum_{\substack{T_j \cap D_j \neq \emptyset \\ T_j \in V_j}} \oint \nabla u \cdot \vec{n}_{ik} \, d\sigma = \\
- \sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} (\phi_i u_i + \phi_j u_j + \phi_k u_k) (n_{ij}^x L_{ij} + n_{ik}^x L_{ik}) \\
- \sum_{\substack{T_j \cap D_j \neq \emptyset, \\ T_j \in V_j}} (\chi_i u_i + \chi_j u_j + \chi_k u_k) (n_{ij}^y L_{ij} + n_{ik}^y L_{ik}),
$$

where $L_{pq} = |Gm_{[pq]}|$ and $\vec{n}_{pq} = (n_{pq}^x \ n_{pq}^y)^\top$ such that $\vec{n}_{pq} \perp \vec{Gm}_{[pq]}$ and $\vec{n}_{pq} \cdot \vec{p}q > 0$. This leads to the following linear system

$$
\begin{align}
\sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} \left(\begin{array}{c} \lambda \\ \mu \end{array}\right) \cdot \left( |Gm_{[ij]}| \vec{n}_{ij} + |Gm_{[ik]}| \vec{n}_{ik} \right) = & \nonumber \\
& - \sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} u_i \left( \phi_i (n_{ij}^x L_{ij} + n_{ik}^x L_{ik}) + \chi_i (n_{ij}^y L_{ij} + n_{ik}^y L_{ik}) \right) \nonumber \\
& - \sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} u_j \left( \phi_j (n_{ij}^x L_{ij} + n_{ik}^x L_{ik}) + \chi_j (n_{ij}^y L_{ij} + n_{ik}^y L_{ik}) \right) \nonumber \\
& - \sum_{\substack{T_i \cap D_i \neq \emptyset, \\ T_i \in V_i}} u_k \left( \phi_k (n_{ij}^x L_{ij} + n_{ik}^x L_{ik}) + \chi_k (n_{ij}^y L_{ij} + n_{ik}^y L_{ik}) \right). \tag{2.5}
\end{align}
$$

Next, we assemble (2.5) for each element $T_i$ results in the global stiffness matrix [S]. Note that each row of this matrix is associated with the vertex $i$ in the mesh and at the $i$th row of [S] it contains contributions of all elements with common vertex $i$.

Next, we formulate the mass matrix in the finite volume discretiza-
tion (2.2). Hence,

$$
\int_{D_i} u \, dx = \sum_{T_i \cap D_i \neq \emptyset} \left( \int_{D_i \cap T_i} u \, dx \right) = \sum_{T_i \cap D_i \neq \emptyset} \left( \int_{T_i \cap V_i} u \, dx + \int_{T_i \cap V_i} u \, dx \right). \quad (2.6)
$$

At the points $m_{[ij]}$ and $G$,

$$
u(m_{[ij]}) = \frac{1}{2}(u_i + u_j), \quad u(G) = \frac{1}{3}(u_i + u_j + u_k).
$$

Thus, the first integral term in (2.6) is approximated as

$$
\begin{align}
\int_{im_{[ij]}G} u \, dx &= \frac{1}{3} |im_{[ij]}G| \left( u_i + \frac{1}{2}(u_i + u_j) + \frac{1}{3}(u_i + u_j + u_k) \right), \nonumber \\
&= \frac{1}{3} |im_{[ij]}G| \left( \frac{11}{6}u_i + \frac{5}{6}u_j + \frac{1}{3}u_k \right). \tag{2.7}
\end{align}
$$

It is straightforward to verify that

$$
\bar{\mathrm{i}} \bar{\mathrm{m}}_{[ij]} = \frac{1}{2} \bar{\mathrm{i}} \bar{\mathrm{j}}, \quad \bar{\mathrm{i}} \bar{\mathrm{G}} = \frac{2}{3} \bar{\mathrm{i}} \bar{\mathrm{m}}_{[jk]}, \quad \bar{\mathrm{i}} \bar{\mathrm{m}}_{[jk]} = \bar{\mathrm{i}} \bar{k} + \frac{1}{2} \bar{\mathrm{k}} \bar{\mathrm{j}}.
$$

Hence,

$$
|im_{[ij]}G| = \frac{1}{2} \left\| \frac{1}{2} \bar{i} \bar{j} \wedge \frac{2}{3} \left( \bar{i} \bar{k} + \frac{1}{2} \bar{j} \bar{k} \right) \right\| = \frac{1}{12} \left\| \bar{i} \bar{j} \wedge \bar{k} \bar{j} \right\| = \frac{1}{6} |T_i|,
$$

and the approximation (2.7) becomes

$$
\int_{im_{[ij]}G} u \, dx = \frac{|T_i|}{18} \left( \frac{11}{6} u_i + \frac{5}{6} u_j + \frac{1}{3} u_k \right).
$$

Similarly, the second integral term in (2.6) is approximated by

$$
\begin{align}
\int_{iGm_{[ik]}} u \, dx &= \frac{1}{3} |iGm_{[ik]}| \left( u_i + \frac{1}{2}(u_i + u_k) + \frac{1}{3}(u_i + u_j + u_k) \right), \nonumber \\
&= \frac{1}{3} |iGm_{[ik]}| \left( \frac{11}{6} u_i + \frac{5}{6} u_j + \frac{1}{3} u_k \right). \tag{2.8}
\end{align}
$$

Again, using the fact that

$$
\bar{\mathrm{i}} \bar{\mathrm{m}}_{[ik]} = \frac{1}{2} \bar{\mathrm{i}} \bar{k}, \\
|iGm_{[ik]}| = \frac{1}{2} \left\| \frac{1}{2} \bar{i} \bar{k} \wedge \frac{2}{3} (\bar{i} \bar{j} + \frac{1}{2} \bar{j} \bar{k}) \right\| = \frac{1}{12} \left\| \bar{i} \bar{k} \wedge \bar{k} \bar{j} \right\| = \frac{1}{6} |T_i|.
$$

the approximation (2.8) becomes

$$
\int_{im_{[ij]}G} u \, dx = \frac{|T_i|}{18} \left( \frac{11}{6} u_i + \frac{1}{3} u_j + \frac{5}{6} u_k \right),
$$

and the approximation of the integral in (2.6) is reduced to

$$
\int_{im_{[ij]}G} u \, dx + \int_{iGm_{[ik]}} u \, dx = \frac{|T_i|}{18} \left( \frac{11}{3} u_i + \frac{7}{6} u_j + \frac{7}{6} u_k \right). \quad (2.9)
$$

The mass matrix [M] can be assembled using the above element matri-
ces in the same manner as performed for the stiffness matrix. It should
be noted that, unlike most of cell-centered finite volume methods where
the assembled mass matrix is diagonal, the mass matrix in the pro-
posed vertex-centered finite volume is dense. In general, cell-centered
methods lead to many more degrees of freedom than vertex-centered
methods on unstructured meshes.

**3. Mesh adaptation procedure**

To increase the efficiency of the proposed vertex-centered finite
volume method we consider a mesh adaptation technique based on a
posteriori error estimates. This class of error estimators can be used
to locally refine the mesh and subsequently minimize the discretiza-
tion errors with an optimal computational cost for wave problems with
sharp gradients and propagating fronts. The error estimates preserve
the main properties of an optimal estimator as upper bound of the com-
puted solution, robustness and local efficiency among others. Similar
error estimates have been widely studied in the literature in both finite
element and finite volume discretizations, see for instance, [2,26,1]. In
the present study, error estimators adapted for the wave problems are
based on a locally postprocessed approximate solution preserving the
conservative flux. Let $u_h$ denote the numerical solution of the system
(2.1) obtained at $t = t_n$, and $S_h$ be a (fine) simplicial mesh obtained from
the associated dual mesh $D_h$, see Fig. 3.1 for an illustration. The approx-
imate flux $-\nabla u_h$ is constant by element and satisfies the local conser-
vation but it does not belong to $\mathbf{H}(\text{div}, \Omega) = \{\mathbf{w} \in L^2(\Omega); \nabla \cdot \mathbf{w} \in L^2(\Omega)\}$
<figure><img src="image_5.png" /><figcaption>Fig. 3.1. Construction of the flux $\theta_h$ by direct prescription.</figcaption></figure>

defined as the space of functions with square-integrable weak divergences. Then, we carry out a flux correction starting from the numerical solution by defining a flux $\theta_h \in H(\text{div}, \Omega)$ which approximates the exact flux and verifies

$$ (\nabla \cdot \theta_h, 1)_{D_i} = (f, 1)_{D_i}, \quad \text{in each control volume } D_i \in D_h, $$

where $(\cdot, \cdot)_{D_i}$ denotes the $L^2$-scalar product on $D_i$. We choose the method of direct prescription to reconstruct $\theta_h$ in the lowest-order Raviart-Thomas-Nédélec space $RTN(S_h)$ which is included in the space $H(\text{div}, \Omega)$, see [26,1] for more details on the reconstruction of $\theta_h$. For instance, on the sub-triangle $iG_1 m_{[i,j]}$, the following system

$$ \theta_h \cdot \vec{n}_1 = -\nabla u_h \cdot \vec{n}_1, $$

$$ \theta_h \cdot \vec{n}_2 = -\frac{1}{2} \left( \nabla u_{h|T_1} \cdot \vec{n}_2 \right) - \frac{1}{2} \left( \nabla u_{h|T_2} \cdot \vec{n}_2 \right), $$

$$ \theta_h \cdot \vec{n}_3 = -\nabla u_h \cdot \vec{n}_3, $$

is solved. Once the flux solution $\theta_h$ is computed, the following a posteriori estimator is established

$$ |||u - u_h||| \le \left( \sum_{D_i \in D_h} (\eta_{R,D_i} + \eta_{DF,D_i})^2 \right)^{\frac{1}{2}}, \quad (3.1) $$

where $||| \cdot |||$ is the energy norm defined as $|||u|||^2 := ||\nabla u||_{\Omega}^2$, $\eta_{DF,D_i}$ the diffusive flux estimator and $\eta_{R,D_i}$ the residual estimator given by

$$ \eta_{DF,D_i} = ||\nabla u_h + \theta_h||_{D_i}^2, \quad \eta_{R,D_i} := m_{D_i} ||f - \nabla \cdot \theta_h||_{D_i}^2, $$

where $m_{D_i}$ is a computable constant depending on the diameter of the dual element $D_i$ and the wave celerity. $m_{D_i}$ is defined as $m_{D_i}^2 = C_p d_{D_i}^2 / c$, with $C_p = 1/\pi^2$. Here, the a posteriori error estimate (3.1) is used as an error indicator along with a given tolerance to adapt the mesh accordingly. For the mesh adaptation, we combine the techniques investigated in [9,5] to develop an efficient adaptive finite volume method for transient wave problems. First, we proceed by refining the mesh using the adaptive algorithm from [5]. Then, the Newest-Vertex-Bisection algorithm proposed in [9] is used for the mesh conformity. It should be noted that there is no more refinement propagation using this combined algorithm compared to the other conventional mesh adaptation techniques as those studied in [5,32]. It should be stressed that the solution is known in the vertex of the primal triangle $T_i$ and the estimator is computed by interpolation in the dual mesh $D_i$ which is obtained by dividing the primal triangle in six sub-triangles.

In the current study, the adaptation procedure is carried out using four steps namely (i) solve the wave problem, (ii) calculate the error estimates, (iii) mark the selected elements, and (iv) refine or/coarsen the current mesh. Hence, given a tolerance $\epsilon$, and let $Estim(T_i)$ denote the normalized error estimate on $T_i$, if $Estim(T_i)$ satisfies the following condition in $T_i$

$$ Estim(T_i) \le \epsilon, $$

then the element $T_i$ is divided into 4 triangles. Multiple refinements can be added by introducing a series of tolerances $\epsilon_m$ with $0 = \epsilon_0 < \epsilon_1 < \dots < \epsilon_m < \epsilon_{m+1} = 1$ and $m \ge 3$. If a macro-element $T_i$ satisfies the condition

$$ \epsilon_m \le Estim(T_i) \le \epsilon_{m+1}, \quad m = 0, 1, \dots, $$

then $T_i$ is divided into $4^m$ triangles. Note that the values of $\{\epsilon_1, \dots, \epsilon_m\}$ can be interpreted as tolerances to be set by the user resulting into a multilevel adaptive procedure. The procedure used in our study for mesh adaptation is based on multilevel refinements and unrefinements allowing the reconstruction of adaptive meshes which dynamically follow the solution of our problem under study. The algorithm begins by calculating the error estimates to make the refinement and unrefinement decisions. Later, we reconstruct a list $\mathcal{L}$ of elements to be refined, their degree of refinement, and those to be unrefined in the computational domain. This step is achieved by filling an integer array denoted for example by $I$ for all triangles of the coarse mesh. At the current time $t = t_n$, and for a macro-element $T_i$, we set $I(T_i) = m$ which means that the element $T_i$ has to be divided into $4^m$ triangles at the next time $t = t_{n+1}$. Hence, starting from a mesh level $r$, made of $N_r$ elements, the next mesh level contains $N_{r+1} = 4 \times N_r$ elements. It is clear that this procedure can be repeated as long as $r < M$ with $M$ is the total number of refinement levels. In order to obtain a mesh with minimum distortion, the procedure opts to divide into two equal parts certain number of edges in the elements to be refined. In Fig. 3.2 we illustrate the procedure for a two-level refinement of triangular elements. The numbers in this figure refer to the number of refinements to be performed for each sub-triangle. Once the division of triangles is done, we proceed to the conformity of the mesh using the Newest-Vertex-Bisection method as follows: for each triangle marked from the primary mesh, a particular edge, called the refinement edge, is selected. The marked element is divided into two new elements by connecting the middle point of the refinement edge to the opposite vertex, which reduces its volume while keeping some regularity of the mesh.

In summary, the adaptive procedure presented in this study starts by computing the error estimates on a fixed coarse mesh. Then, marking the nodes concerned by the refinement process using the provided tolerances. Next, the interval is split into subintervals which give the decision on how many times each cell of the reference mesh should be divided. This procedure satisfies the condition that the current level of refinement for the neighboring cells cannot be higher than one. At this stage, the conformity is ensured by the Newest-Vertex-Bisection algorithm as illustrated in Fig. 3.2. After a given number of iterations, the new criterion is computed on the reference mesh and a new refined mesh is created. An interpolation procedure is used to transfer the solutions from the old mesh to the new refined mesh. These algorithmic steps are summarized in Algorithm 1.
<figure><img src="image_6.png" /><figcaption>Fig. 3.2. Illustration of a two-level refining for triangular elements. The numbers in the figure refer to the number of refinements to be performed for each sub-triangle.</figcaption></figure>

**Algorithm 1** Adaptive procedure used in the present study.

Require: $\mathcal{N}_{Ref}$ maximum level of refinements

Require: $\mathcal{N}_{Lev}$ number of multilevel refinements

1: **for** $k = 1 : \mathcal{N}_{Ref}$ **do**

2: Create the mesh data structures

3: Compute the numerical solution

4: Calculate the error estimator

5: Set the criteria of adaptivity

6: **for** $l = 1 : \mathcal{N}_{Lev}$ **do**

7: Preform the mesh refinement using the adaptive procedure

8: Ensure the conformity using the Newest-Vertex-Bisection algorithm

9: Mark the new triangles created to be refined

10: **end for**

11: **end for**

## 4. Time integration

The semi-discretized equations of an initial-value problem of the form (2.2) can be rewritten in a common notation as

$$[\mathbf{M}] \{\ddot{\mathbf{U}}\} + [\mathbf{C}] \{\dot{\mathbf{U}}\} + [\mathbf{K}] \{\mathbf{U}\} = \{\mathbf{F}\}, \quad (4.1)$$

where $[\mathbf{M}]$ represents the mass matrix, $[\mathbf{C}]$ the damping matrix, $[\mathbf{K}]$ the stiffness matrix, $\{\mathbf{U}\}$ the nodal solution vector and $\{\mathbf{F}\}$ the external force vector at time $t$. A superposed dot over a quantity denotes its derivative with respect to time. Numerical time integration of the system (4.1) is often carried out using a Newmark time stepping scheme [36]. To discrete the time domain, the time interval is divided into small subintervals $[t_n, t_{n+1}]$ with stepsize $\Delta t = t_{n+1} - t_n$ for $n = 0, 1, \dots$. We use the notation $\mathbf{W}^n$ to denote the value of the solution $\mathbf{W}(t)$ at time $t = t_n$. The Newmark scheme is based on the following approximations

$$\begin{aligned} \{\mathbf{U}^{n+1}\} &= \{\mathbf{U}^n\} + \Delta t \{\dot{\mathbf{U}}^n\} + \frac{\Delta t^2}{2} \left( (1-2\beta) \{\ddot{\mathbf{U}}^n\} + 2\beta \{\ddot{U}^{n+1}\} \right), & (4.2) \\ \{\ddot{\mathbf{U}}^{n+1}\} &= \{\dot{\mathbf{U}}^n\} + (1-\gamma) \Delta t \{\dot{\mathbf{U}}^n\} + \gamma \Delta t \{\ddot{U}^{n+1}\}, \end{aligned}$$

where $\gamma$ and $\beta$ are the Newmark parameters which determine the stability and accuracy characteristics of the integration scheme under consideration. The Newmark method is unconditionally stable when $\gamma$ is greater than $\frac{1}{2}$ and $\beta$ is greater than $(\frac{1}{2} + \gamma)^2/4$. The most widely used choice is $\gamma = \frac{1}{2}$ and $\beta = \frac{1}{4}$, which is unconditionally stable and it does not introduce any artificial damping into the solution. This choice is also adopted in the present study. For the solution $\mathbf{U}^{n+1}$ at time $t = t_n + \Delta t$, the equilibrium system is recurrently solved as

$$[\mathbf{K}^*] \{\mathbf{U}^{n+1}\} = \{\mathbf{F}^*\}, \quad (4.3)$$

where the effective stiffness matrix $[\mathbf{K}^*]$ and effective force vector $\{\mathbf{F}^*\}$ are defined by

$$\begin{aligned} [\mathbf{K}^*] &= \left[ \frac{1}{\beta \Delta t^2} [\mathbf{M}] + \frac{\gamma}{\beta \Delta t} [\mathbf{C}] + [\mathbf{K}] \right], \\ \{\mathbf{F}^*\} &= \{\mathbf{F}^n\} + \left[ \frac{1}{\beta \Delta t^2} [\mathbf{M}] + \frac{\gamma}{\beta \Delta t} [\mathbf{C}] \right] \{\mathbf{U}^n\} \end{aligned}$$

$$+ \left[ \frac{1}{\beta \Delta t} [\mathbf{M}] + \left( \frac{\gamma}{\beta} - 1 \right) [\mathbf{C}] \right] \{\ddot{\mathbf{U}}^n\} + \left[ \frac{1 - 2\beta}{2\beta} [\mathbf{M}] - \Delta t \left( 1 - \frac{\gamma}{2\beta} \right) [\mathbf{C}] \right] \{\dot{\mathbf{U}}^n\}.$$

In summary, Algorithm 2 presents the steps used by the Newmark scheme for the time integration. Note that most of the computational effort is required by the solution of the linear system (4.3). Here, the system matrix may be decomposed at each adaptive step and retained to be reused with just updating the right-hand side of the linear system of equations at all the following time steps within the fixed mesh. In our simulations on fixed meshes, the matrix is decomposed into an $\mathbf{L}\mathbf{U}\mathbf{L}^\top$ factorization, then the solution is reduced to backward/forward substitutions after updating the right-hand side vector at every time step.

**Algorithm 2** Newmark scheme for the time integration.

1: Assemble the finite volume matrices $[\mathbf{M}]$ and $[\mathbf{K}]$.

2: Define the damping matrix $[\mathbf{C}]$.

3: **for** $n = 0$ **do**

4: The initial solutions $\{\mathbf{U}^0\}$ and $\{\ddot{\mathbf{U}}^0\}$ are given. Calculate the initial acceleration as

$$\{\ddot{\mathbf{U}}^0\} = [\mathbf{M}]^{-1} \{[\mathbf{F}] - [\mathbf{C}] \{\dot{\mathbf{U}}^0\} - [\mathbf{K}] \{\mathbf{U}^0\}\}.$$

5: Calculate the integration constants

$$b_1 = \frac{1}{\beta \Delta t^2},$$

$$b_2 = -\frac{1}{\beta \Delta t},$$

$$b_3 = \frac{2\beta - 1}{2\beta},$$

$$b_4 = \gamma \Delta t b_1,$$

$$b_5 = 1 + \gamma \Delta t b_2,$$

$$b_6 = \Delta t (1 + \gamma b_3 - \gamma).$$

6: Form the constant matrices

$$[\mathbf{K}^*] = b_1 [\mathbf{M}] + b_4 [\mathbf{C}] + [\mathbf{K}], \quad [\mathbf{A}_1] = b_1 [\mathbf{M}] + b_4 [\mathbf{C}],$$

$$[\mathbf{A}_2] = b_2 [\mathbf{M}] + b_5 [\mathbf{C}], \quad [\mathbf{A}_3] = b_3 [\mathbf{M}] + b_6 [\mathbf{C}].$$

7: Calculate the inverse matrix $[\mathbf{K}^*]^{-1}$.

8: **end for**

9: **for** $n = 0, 1, 2, \dots$ **do**

10: Calculate the effective force

$$\{[\mathbf{F}^*]\} = \{[\mathbf{F}^n]\} + [\mathbf{A}_1] \{[\mathbf{U}^n]\} - [\mathbf{A}_2] \{\ddot{\mathbf{U}}^n\} - [\mathbf{A}_3] \{\dot{\mathbf{U}}^n\}.$$

11: Update the new solution

$$\{\ddot{\mathbf{U}}^{n+1}\} = [\mathbf{K}^*]^{-1} \{\mathbf{F}^*\}.$$

12: Update the new velocity and acceleration solutions as

$$\{\ddot{\mathbf{U}}^{n+1}\} = b_4 \{[\mathbf{U}^{n+1}] - \{\mathbf{U}^n\}\} + b_5 \{\ddot{\mathbf{U}}^n\} + b_6 \{\dot{\mathbf{U}}^n\},$$

$$\{\dot{\mathbf{U}}^{n+1}\} = b_1 \{[\mathbf{U}^{n+1}] - \{\mathbf{U}^n\}\} + b_2 \{\mathbf{U}^n\} + b_3 \{\ddot{\mathbf{U}}^n\}.$$

13: **end for**
<table><caption>**Table 5.1**
Results for $L^1$-error and $L^2$-error obtained for the propagation of a radial wave using FEM and FVM at the time $t = 0.5$ on different meshes.</caption><thead><tr><th rowspan="3">Mesh</th><th rowspan="3"># ele</th><th rowspan="3"># node</th><th colspan="4">FEM</th><th colspan="4">k = 4π</th></tr><tr><th colspan="2">k = 2π</th><th colspan="2">k = 4π</th><th colspan="2">k = 2π</th><th colspan="2">k = 4π</th></tr><tr><th>L<sup>1</sup>-error</th><th>q</th><th>L<sup>2</sup>-error</th><th>q</th><th>L<sup>1</sup>-error</th><th>q</th><th>L<sup>2</sup>-error</th><th>q</th></tr></thead><tbody><tr><td>h<sub>1</sub></td><td>787</td><td>431</td><td>1.28E-02</td><td>—</td><td>1.40E-02</td><td>—</td><td>6.22E-02</td><td>—</td><td>6.73E-02</td><td>—</td></tr><tr><td>h<sub>1/2</sub></td><td>1541</td><td>830</td><td>6.59E-03</td><td>0.958</td><td>7.14E-03</td><td>0.971</td><td>3.52E-02</td><td>0.821</td><td>3.73E-02</td><td>0.851</td></tr><tr><td>h<sub>1/4</sub></td><td>3134</td><td>1639</td><td>3.30E-03</td><td>0.997</td><td>3.57E-03</td><td>0.999</td><td>1.60E-02</td><td>1.137</td><td>1.73E-02</td><td>1.108</td></tr><tr><td>h<sub>1/8</sub></td><td>6269</td><td>3259</td><td>1.71E-03</td><td>0.948</td><td>1.84E-03</td><td>0.956</td><td>7.94E-03</td><td>1.011</td><td>8.64E-03</td><td>1.002</td></tr></tbody></table><table><thead><tr><th rowspan="3">Mesh</th><th rowspan="3"># ele</th><th rowspan="3"># node</th><th colspan="4">FVM</th><th colspan="4">k = 4π</th></tr><tr><th colspan="2">k = 2π</th><th colspan="2">k = 4π</th><th colspan="2">k = 2π</th><th colspan="2">k = 4π</th></tr><tr><th>L<sup>1</sup>-error</th><th>q</th><th>L<sup>2</sup>-error</th><th>q</th><th>L<sup>1</sup>-error</th><th>q</th><th>L<sup>2</sup>-error</th><th>q</th></tr></thead><tbody><tr><td>h<sub>1</sub></td><td>787</td><td>431</td><td>1.13E-02</td><td>—</td><td>1.36E-02</td><td>—</td><td>5.72E-02</td><td>—</td><td>6.49E-02</td><td>—</td></tr><tr><td>h<sub>1/2</sub></td><td>1541</td><td>830</td><td>4.86E-03</td><td>1.217</td><td>5.74E-03</td><td>1.245</td><td>2.28E-02</td><td>1.324</td><td>2.60E-02</td><td>1.317</td></tr><tr><td>h<sub>1/4</sub></td><td>3134</td><td>1639</td><td>1.85E-03</td><td>1.391</td><td>2.20E-03</td><td>1.380</td><td>8.21E-03</td><td>1.476</td><td>9.32E-03</td><td>1.482</td></tr><tr><td>h<sub>1/8</sub></td><td>6269</td><td>3259</td><td>6.70E-04</td><td>1.467</td><td>8.05E-04</td><td>1.454</td><td>2.20E-03</td><td>1.896</td><td>2.54E-03</td><td>1.875</td></tr></tbody></table>

## 5. Results and examples

A number of numerical examples are selected to illustrate the accuracy of the new adaptive vertex-centered unstructured finite volume method introduced in the above sections. In the first example, we focus on the advantages of the method compared to the standard finite element method using fixed meshes. Then, we examine the performance of the proposed adaptive approach in a problem of wave propagation in a domain without and with embedded barriers where the waves scatter, diffract and reflect. For the evaluation of adaptive results, reference solutions are generated using very fine fixed meshes. In all the simulations presented in this section, we used $c = 1$ and $h = 1$ in the wave problem (2.1) and fixed time steps $\Delta t$. In principle, no limits are required for the time step $\Delta t$ as the Newmark scheme is implicit by reconstruction and unconditionally stable. Therefore, the selection of time steps in our simulations is mainly carried out based on the accuracy purposes and not on stability restrictions. It should be noted that although the FVM and the FEM can have similar flexibilities in dealing with complex geometries and/or heterogeneous domains, the FVM remains to be a less popular choice for solving wave problems. Only limited comparisons are available in the literature between the two methods for the solution of the wave problem (2.1). The second example is used for the numerical assessment of the mesh adaptation used in FVM for solving wave propagations. All the computations were performed on an Intel® Core(TM) i7-7700HQ CPU @ 2.80 GHz with 8 GB of RAM.

### 5.1. Propagation of a radial wave

In the first test case we aim to recover a circular wave that is propagating from the origin of the coordinate system towards the infinity as studied in [17]. The exact solution of this wave problem is defined by

$$u_{\text{exact}}(t, x, y) = \exp(kr - \omega t),$$

where $r$ is the radial coordinate defined by $r = \sqrt{x^2 + y^2}$. We aim to recover the wave inside a unit squared computational domain $\Omega = [1, 2] \times [1, 2]$. To avoid any errors resulting from artificial boundary conditions we impose the analytical solution on the boundary $\Gamma$ using the conditions (2.1b). The initial conditions (2.1c) and (2.1d) are also evaluated using the analytical solution. The angular frequency in this example is fixed at $\omega = 2\pi$ while two different wavenumbers, namely, $k = 2\pi$ and $k = 4\pi$ are considered. The relative error norms $L^1$-error and $L^2$-error for the numerical solution are evaluated at time $t_n$ as

$$L^1\text{-error} = \frac{\int_{\Omega} |u_h^n - u_{\text{exact}}^n| d\mathbf{x}}{\int_{\Omega} |u_{\text{exact}}^n| d\mathbf{x}}, \quad L^2\text{-error} = \frac{\left( \int_{\Omega} |u_h^n - u_{\text{exact}}^n|^2 d\mathbf{x} \right)^{\frac{1}{2}}}{\left( \int_{\Omega} |u_{\text{exact}}^n|^2 d\mathbf{x} \right)^{\frac{1}{2}}},$$

where $u_{\text{exact}}^n$ and $u_h^n$ are, respectively, the exact and numerical solutions at a gridpoint $x_h$ and the time $t_n$.

First the domain is discretised with an unstructured mesh into 787 linear triangular elements and 431 nodes. Then three consecutive mesh refinements are applied. The resulting four meshes are referred to as $h_1$ for the first mesh and $h_1, h_1, h_1$ for the refinements. The total number of elements and nodes are listed in Table 5.1 for individual meshes. The four meshes are used to solve the problem for the two considered wavenumbers using the FEM. The same meshes are used to solve the problem again but using the FVM. The considered timestep size in all the computations is $\Delta t = 0.005$. Table 5.1 summarizes the $L^2$-error and $L^1$-error for both FEM and FVM at the instant $t = 0.5$. Table 5.1 also includes the convergence rate $q$ for each method obtained for the considered mesh refinements. The results in this table show that the FVM consistently outperform the FEM in terms of errors and convergence rates. In all the presented results using the same mesh, the FVM always leads to smaller errors. Note that since linear elements are considered in these simulations the expected FEM convergence rate is 1. The results show that the actual FEM convergence rate is slightly smaller than 1 in most cases which could be attributed to the nature of the manufactured solutions and the values of wavenumbers used. However, the FVM convergence rate is consistently larger than 1 and in the case of $k = 4\pi$ and $h_1$ the convergence rate is $q \approx 1.9$. The results in this section show the potential saving in the number of degrees of freedom that can be achieved when using the proposed FVM in solving wave problems. This is particularly important in this class of problems as in general wave problems require a large number of degrees of freedom especially at high wave numbers.

### 5.2. Wave reflections in a rectangular domain

To study the performance of the proposed adaptive approach we next consider a vertical Gaussian pulse in a rectangular domain. The pulse will result in two vertical Gaussian plane waves moving from the pulse center towards the domain ends without creating any further oscillations. The plane waves then hit the domain sides and reflect inward to meet again. Note that, for the adaptive approach to be successful it needs to refine the mesh only where sharp gradients at high-amplitude
<figure><img src="image_7.png" /><figcaption>Fig. 5.1. Reference solutions obtained for the wave reflections in a rectangular domain using the FVM on the reference mesh with 238555 elements and 120037 nodes at three different times $t = 2, 5$ and $7$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5.2. Results obtained for the wave reflections in a rectangular domain using FEM (first row) and FVM (second row) on a fixed coarse mesh with 15761 nodes at three different times $t = 2, 5$ and $7$.</figcaption></figure>

waves are detected in the computational domain at the current time, otherwise coarse elements are enough to resolve the constant parts of the solution. These wave features are very interesting to check the accuracy of the criteria (3.1) considered in this study for mesh adaptation in the finite volume solution of time-dependent wave reflections. Hence, the aim in this example is to validate the ability of the proposed adaptive approach to accurately capture and closely follow the correct wave locations in the domain. The problem is also solved using the linear finite elements to be compared to the proposed FVM on fixed meshes. The FVM is then used together with the proposed mesh adaptation procedure and comparisons to results obtained using fixed meshes are presented. Here, the considered computational domain is defined as a rectangle $\Omega = [0, 2] \times [0, 1]$ and a damping term is added to the dynamic system (4.1) as

$$[C] = \sigma [M].$$

In our simulations, the wave speed $c$ and the damping coefficient $\sigma$ are assumed to be constants as $c = 1$ and $\sigma = 0.1$. The pulse is centered at $x_0 = 0.8$, its width is $\delta = 0.05$ and reflective conditions are imposed on the domain boundary $\Gamma$, so that a wave is fully reflected if it hits the boundary. The initial conditions producing the pulse are given as

$$u_0(x, y) = \exp\left(-\frac{(x-x_0)^2}{\delta^2}\right), \quad v_0(x, y) = 0. \qquad (5.1)$$

First, a coarse mesh of 31024 linear elements and a total of 15761 nodes, is used in both the FEM and the FVM. The first mesh is referred to as Mesh 1. This is then repeated with Mesh 2 composed of 62073 elements and 31419 nodes. Next, the adaptive FVM is used starting from a coarse initial mesh of 4438 elements and 2316 nodes. We consider single-level, two-level and three-level adaptive procedure using the tolerances $\epsilon_1 = 0.25$, $\epsilon_2 = 0.5$ and $\epsilon_3 = 0.75$. To validate these solutions a reference solution is also computed using a highly refined fixed mesh of

238555 elements and 120037 nodes. In all these cases, a small time step $\Delta t = 5 \times 10^{-3}$ is considered and the simulations are carried out for the final time $t = 7$. It should be noted that the time step $\Delta t$ is chosen so small such that the time integration errors are minimized in the reference solution level. For the consistency in the compared results, the same time step $\Delta t$ is also used in the other considered meshes. Fig. 5.1 exhibits the snapshots of the reference solutions obtained at three different instants $t = 2, t = 5$ and $t = 7$. As expected the pulse initiated at $t = 0$ in $x = 0.8$ splits into two wavefronts, one moving to the left and the other to the right end of the domain with the speed $c = 1$. The wavefronts hit the domain vertical sides and reflect back to merge into a single wave at $t = 2$ in $x = 1.2$ which can be clearly seen in the results shown in Fig. 5.1. The wavefronts pass each other moving in opposite directions to hit the domain ends and reflected again to meet at the initial pulse location i.e. in $x = 0.8$ at $t = 4$. This pattern is then repeated every 4 time units until the pulse is damped away in the computational domain. It is clear that the dynamics of the solution is correctly recovered using the FVM on the reference mesh.

Fig. 5.2 compares solutions obtained using the FEM to those obtained using the FVM both on the fixed Mesh 1. It is clear that both methods seem to recover the problem dynamics on the considered mesh. However, the FEM seems to also recover a shadow of non-physical oscillations following each wavefront in the computational domain. The shadow becomes longer and more pronounced as the simulation continues which it can be seen when comparing the solution snapshots at $t = 5$ and $t = 7$ in Fig. 5.2. Furthermore, at $t = 7$ the wavefront seems to be contaminated with numerical errors so that the wavefront does not remain straight as expected. When compared to the FEM solution, we can see that the FVM solution follows a similar pattern. However, the FVM seems to lead to a higher accuracy as the shadow region seems to be shorter than in the FEM solutions. In addition, the FVM solution at $t = 7$ seems to be more stable where the wavefront remains visually straight in the computational domain. In
<figure><img src="image_9.png" /><figcaption>Fig. 5.3. Results obtained for the wave reflections in a rectangular domain using FEM (first row) and FVM (second row) on a fixed fine mesh with 31419 nodes at three different times $t = 2, 5$ and $7$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5.4. Adaptive meshes obtained for FVM solving the wave reflections in a rectangular domain using single-level adaptation (first row), two-level adaptation (second row) and three-level adaptation (third row) at three different times $t = 2, 5$ and $7$.</figcaption></figure>

Fig. 5.3 we repeat the same comparisons between both methods but on the Mesh 2 formed with about double nodes and elements of Mesh 1. The same observations made on the previous figure are still correct in general for the computed solutions. However, and as expected, the error seems to be reduced using both FEM and FVM on the considered fixed meshes. Again one may observe that the errors computed using the FVM remain smaller compared to those computed using the FEM.

To evaluate the performance of the proposed adaptive FVM for this wave problem, we illustrate in Fig. 5.4 the adaptive meshes obtained with each adaptive level at the same instants as those displayed in Fig. 5.1. For a better insight, the refinement produced by each adaptive level is shown with a different color. Notice that all adaptive levels are initialized starting from the same mesh composed of 4438 elements and 2316 nodes. The mesh is then refined at every timestep based on

the proposed adaptive algorithm according to single-level, two-level and three-level adaptations using the selected tolerances. The number of nodes and elements in each of the considered three levels and at the considered instants $t = 2, 5$ and $7$ are summarized in Table 5.2. It should be also noted that, due to the criterion used for the mesh adaptivity in this study, no regularization processes as those widely used in finite element methods are needed. Fig. 5.4 demonstrates that the mesh is only refined at the right location of the Gaussian waves and coarse elements are used elsewhere. As the waves propagate, the refined parts of the mesh travel as well and closely follow the location of the wavefronts. Moreover, when moving from the single-level to two-level adaptations, the mesh is refined further in the area close to the center of the wave. This is a reasonable behavior as the considered waves have steeper gradients closer to the center. Indeed this behavior is also observed when
<caption>**Table 5.2**</caption>

Mesh statistics obtained for the adaptive FVM solving the wave reflections in a rectangular domain at three different times using single-level, two-level and three-level procedures. The CPU times are given in seconds.

<table><thead><tr><th rowspan="2"></th><th colspan="3">t = 2</th><th colspan="3">t = 5</th><th colspan="3">t = 7</th></tr><tr><th># nodes</th><th># elements</th><th>CPU</th><th># nodes</th><th># elements</th><th>CPU</th><th># nodes</th><th># elements</th><th>CPU</th></tr></thead><tbody><tr><td>Single-level</td><td>4159</td><td>8093</td><td>8.80</td><td>5705</td><td>11160</td><td>28.70</td><td>5677</td><td>11101</td><td>40.03</td></tr><tr><td>Two-level</td><td>7253</td><td>14258</td><td>15.71</td><td>14401</td><td>28471</td><td>77.75</td><td>12321</td><td>24341</td><td>107.72</td></tr><tr><td>Three-level</td><td>9911</td><td>19546</td><td>21.15</td><td>22212</td><td>42065</td><td>121.44</td><td>22663</td><td>44973</td><td>179.81</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 5.5.** Cross-sections of the numerical solutions for the wave reflections in a rectangular domain at $y = 0.5$ obtained using FEM on fixed meshes (first row) and FVM on fixed and adaptive meshes (second row) at three different times $t = 2, 5$ and $7$.</figcaption></figure>

switching from two-level to three-level adaptations where the mesh is only refined closer to the center line. Note that one could consider not to run the adaptation for each time step but in longer time windows to reduce the computational cost. This is possible in our implementation but we decided to run it for each time step for two reasons: (i) to demonstrate the potential of the proposed techniques, and (ii) to avoid missing wave patterns which may occur in those longer time windows in the adaptivity.

To have a better insight into the problem, we show in Fig. 5.5 cross-sections of the computed solutions at $y = 0.5$ and at the same instants considered before i.e. $t = 2, 5$ and $7$. Here, we also present the results obtained by the FEM along side with those obtained using the proposed FVM on fixed meshes and those obtained using FVM on adaptive meshes using single-level, two-level and three-level adaptive procedures. Note that both sets of results are plotted against the reference solution. The FEM results exhibit a clear deviation from the reference solution. It is also clear from these cross-section plots that the results on the fixed Mesh 1 are significantly less accurate than the reference solutions when compared to those results obtained on the finer fixed Mesh 2. Nevertheless, both meshes generate non-physical oscillations in the wave shadow regions and underestimate the wave amplitudes. This behavior becomes more pronounced when the solution progresses in time. Similarly the FVM on fixed meshes seems to lead to spurious oscillations and amplitude underestimation but the oscillations are of a smaller magnitude and the underestimation is also less significant than those in the FEM re-

sults. However, this behavior is largely reduced when considering FVM on adapted meshes using the single-level and two-level adaptive procedures. It is also worth mentioning that no stabilization techniques have been used in the FEM. In addition, considering the three-level adaptive procedure, it can be seen that the computed results practically match those of the reference solutions. These results are possible despite the significant reduction in the total number of degrees of freedom from about 120000 in the reference solution to around 23000 in the results obtained using the three-level adaptive procedure. Indeed the number of degrees of freedom can be much less in the adaptive approach at $t = 2$ as the two wavefronts merge together into one wave and the number of degrees of freedom reduces to less than 10000 as can be seen in Table 5.2.

### 5.3. Wave scattering and diffraction

Our last example is the same as the previous example but two thin rectangular barriers forming a narrow gap are included in the computational domain so that the wave is scattered by the plates and is diffracted around the gap. As for the external boundaries, homogeneous Neumann conditions are also imposed on the internal barriers. It should be noted that these obstacles represent a major challenge for numerical methods such as the FEM and FVM due to the singularities appeared at the barriers corners. Furthermore, capturing detailed structures of the wave scattering and diffraction at the narrow gap requires a
<figure><img src="image_12.png" /><figcaption>Fig. 5.6. Coarse fixed mesh (left), fine fixed mesh (middle) and initial adaptive mesh (right) used for the wave scattering and diffraction.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 5.7. Results obtained for the wave scattering and diffraction on the coarse mesh (first column), on the fine mesh (second column) and on the adaptive mesh (third column) at four different times. From top to bottom $t = 0.2, 0.4, 0.65$ and 1.</figcaption></figure>

highly refined mesh. Similar to the previous problem, the same computational domain is considered again $\Omega = [0, 2] \times [0, 1]$. The wave speed, the damping coefficient and the time step are also unchanged i.e. $c = 1$, $\sigma = 0.1$ and $\Delta t = 10^{-3}$. However, we now consider a shorter wave width $\delta = 0.009$. It should be noted that a similar problem has been solved in [24] using a Runge-Kutta time stepping scheme.

In our simulations for this test example we only consider the adaptive FVM using three-level adaptive procedure using the tolerances $\epsilon_1 = 0.25$, $\epsilon_2 = 0.5$ and $\epsilon_3 = 0.75$. The initial mesh to be used in the adaptation is selected to be coarse composed of 4546 nodes and 8732 elements. For a comparison reasons, we also present numerical results obtained using the FVM on two fixed meshes. The first mesh is relatively coarse but much finer than the previous one used as initial mesh for the adaptation and it is composed of 15604 nodes and 30552 elements while, the second fixed mesh is highly refined and composed of

31017 nodes and 61027 elements, see Fig. 5.6 for an illustration of these meshes. The obtained solutions using the considered fixed and adaptive meshes are shown in Fig. 5.7 for four different instants $t = 0.2, 0.4, 0.65$ and 1. The obtained adaptive meshes at the corresponding instants are illustrated in Fig. 5.8. As it can be seen from the wave solutions in Fig. 5.7, the results obtained using the adaptive FVM with three-level adaptive procedure are significantly more accurate than those obtained using the FVM on fixed meshes. The spurious oscillations detected in the solution obtained on the coarse fixed mesh and to an extent on the fine fixed mesh completely disappear in those results obtained using adaptive meshes. As expected, the two generated waves propagate in the computational domain generating complex patterns once reaching the two obstacles. These wave patterns are very well captured by the proposed adaptive finite volume method and the small wave features at the gap between the two barriers are also well resolved using
<figure><img src="image_14.png" /><figcaption>Fig. 5.8. Adaptive meshes obtained for the wave scattering and diffraction at the four considered times. From top to bottom $t = 0.2, 0.4, 0.65$ and 1.</figcaption></figure>

<caption>**Table 5.3**</caption>

Mesh statistics and computational costs obtained for the adaptive FVM solving the wave scattering and diffraction at four different times the three-level procedure.

<table><thead><tr><th></th><th># nodes</th><th># elements</th><th>CPU time (in minutes)</th></tr></thead><tbody><tr><td>t = 0.2</td><td>29749</td><td>15143</td><td>1.732</td></tr><tr><td>t = 0.4</td><td>32345</td><td>16464</td><td>3.276</td></tr><tr><td>t = 0.65</td><td>43941</td><td>22294</td><td>4.712</td></tr><tr><td>t = 1</td><td>53938</td><td>27290</td><td>6.508</td></tr></tbody></table>

our method. Here, the initial mesh is intentionally selected to be sufficiently coarse, being composed of 31024 elements only. At early times of the simulation, the wave dispersion zone is small so that elements with large aspect ratio are not generated in the adapted meshes. Notice that small elements are needed in the gap area in order to capture the diffraction of the wave and its scattering. These results also exhibit the scattered waves formed by the reflection against the obstacles move back in the wave direction. The scattering of the waves around the obstacles and the gap zone can also be identified in the adapted meshes in Fig. 5.8. It is also evident that the considered adaptive mesh procedure is capable to capture the complex features of the wave with a high level of accuracy. For instance, the comparison with results reported in [24] for a similar problem shows that the proposed adaptive finite volume method is able to predict complex wave interactions with high accuracy.

To quantify the results for this test problem we summarize in Table 5.3 the mesh statistics and the CPU times for the considered adaptive meshes using the three-level adaptive procedure at the four considered times $t = 0.2, 0.4, 0.65$ and 1. As would be expected, the denser the mesh the greater the accuracy of the numerical predictions and the computational effort required. Furthermore, the simulations carried out on the adaptive meshes are of the same order of accuracy as those on the fine mesh despite the fact that the adaptive mesh uses significantly fewer degrees of freedom. For the considered wave conditions, using the error estimate (3.1) in the adaptive criteria seems to guarantee a high accuracy and efficiency in the computed results. Computational results obtained for the considered wave problems reveal that simulations using the proposed adaptive finite volume method are more accurate than those using fine fixed meshes. In addition, the developed adaptation procedure gives encouraging results, with global errors always smaller than those obtained by uniformly refining the whole mesh.

## 6. Conclusions

We proposed a new vertex-centered unstructured finite volume method for solving two-dimensional time-domain wave problems on adaptive meshes. The new method incorporates the combination of a vertex-centered finite volume for the space discretization and an implicit Newmark scheme for the time integration. As a criterion for the mesh adaptation, we used an error estimator in the energy norm and multilevel adaptive procedures have been implemented. The proposed vertex-centered finite volume method is suitable for complex geometries, independent of the sizes and arrangement of the mesh elements, and it can be performed using time steps larger than those required for its explicit counterparts. In addition, the proposed method is simple, stable and eliminates the numerical difficulties related to the fine spatial discretization required to capture the solution structures in wave problems. Numerical results for examples of wave propagation, reflection, scattering and diffraction show that the method provides very accurate results. A comparison to the conventional finite element method and other reference solutions shows the capabilities of the current vertex-centered finite volume method to solve the time-dependent wave problems. Moreover, the presented results demonstrate the capability of the adaptive vertex-centered finite volume method to provide insight to complex wave features. It has been shown that the proposed adaptive vertex-centered unstructured finite volume method enjoys the computational advantages and achieves accurate solutions for wave propagations. For a future work we aim to extend the work presented here to solving the transient wave equations in three space dimensions using unstructured tetrahedral meshes. Moreover, to further improve the efficiency of the proposed method, it will also be of interest to develop space-time adaptivity and investigate enriching the approximation space with oscillatory basis functions.

## Data availability

Data will be made available on request.

## References

[1] M. Ainsworth, T. Vejchodský, Robust error bounds for finite element approximation of reaction-diffusion problems with non-constant reaction coefficient in arbitrary space dimension, Comput. Methods Appl. Mech. Eng. 281 (2014) 184–199.
[2] M. Ainsworth, T. Vejchodský, A simple approach to reliable and robust a posteriori error estimation for singularly perturbed problems, Comput. Methods Appl. Mech. Eng. 353 (2019) 373–390.

[3] I. Babuška, S.A. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, SIAM J. Numer. Anal. 34 (6) (1997) 2392–2423.

[4] Jeffrey W. Banks, William D. Henshaw, Upwind schemes for the wave equation in second-order form, J. Comput. Phys. 231 (17) (2012) 5854–5889.

[5] F. Benkhaldoun, I. Elmahi, M. Seaid, Well-balanced finite volume schemes for pollu- tant transport by shallow water equations on unstructured meshes, J. Comput. Phys. 226 (1) (2007) 180–203.

[6] S. Bilbao, B. Hamilton, Higher-order accurate two-step finite difference schemes for the many-dimensional wave equation, J. Comput. Phys. 367 (2018) 134–165.

[7] D. Borisov, R. Modrak, F. Gao, J. Tromp, 3D elastic full-waveform inversion of sur- face waves in the presence of irregular topography using an envelope-based misfit function, Geophysics 83 (1) (2018) R1–R11.

[8] D. Braess, J. Schöberl, Equilibrated residual error estimator for edge elements, Math. Comput. 77 (262) (2008) 651–672.

[9] L. Chen, C. Zhang, A coarsening algorithm on adaptive grids by newest vertex bi- section and its applications, J. Comput. Math. (2010) 767–789.

[10] L.L. Chen, Y. Zhang, H. Lian, E. Atroshchenko, C. Ding, S.P.A. Bordas, Seamless integration of computer-aided geometric modeling and acoustic simulation: isoge- metric boundary element methods based on Catmull-Clark subdivision surfaces, Adv. Eng. Softw. 149 (2020) 102879.

[11] Z. Chen, F. Zhang, J. Liu, B. Chen, A vertex-based reconstruction for cell-centered finite-volume discretization on unstructured grids, J. Comput. Phys. (2021) 110827.

[12] A. Deraemaeker, I. Babuška, P. Bouillard, Dispersion and pollution of the FEM so- lution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Methods Eng. 46 (4) (1999) 471–499.

[13] C.B. Dilgen, N. Aage, Generalized shape optimization of transient vibroacous- tic problems using cut elements, Int. J. Numer. Methods Eng. 122 (6) (2021) 1578–1601.

[14] G.C. Diwan, M.S. Mohamed, Pollution studies for high order isogeometric analysis and finite element for acoustic problems, Comput. Methods Appl. Mech. Eng. 350 (2019) 701–718.

[15] G.C. Diwan, M.S. Mohamed, Iterative solution of Helmholtz problem with high- order isogeometric analysis and finite element method at mid-range frequencies, Comput. Methods Appl. Mech. Eng. 363 (2020) 112855.

[16] G.C. Diwan, M.S. Mohamed, Iterative solution with shifted Laplace preconditioner for plane wave enriched isogeometric analysis and finite element discretization for high-frequency acoustics, Comput. Methods Appl. Mech. Eng. 384 (2021) 114006.

[17] M. Drolia, M.S. Mohamed, O. Laghouche, M. Seaid, A. El Kacimi, Explicit time integration with lumped mass matrix for enriched finite elements solution of time domain wave problems, Appl. Math. Model. 77 (2020) 1273–1293.

[18] M. Drolia, M.S. Mohamed, O. Laghouche, M. Seaid, J. Trevelyan, Enriched finite elements for initial-value problem of transverse electromagnetic waves in time do- main, Comput. Struct. 182 (2017) 354–367.

[19] M. Dumbser, M. Käser, Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems, J. Comput. Phys. 221 (2) (2007) 693–723.

[20] M. Dumbser, M. Käser, J. De La Puente, Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D, Geophys. J. Int. 171 (2) (2007) 665–694.

[21] M. Duprez, S.P.A. Bordas, M. Bucki, H.P. Bui, F. Chouly, V. Lleras, C. Lobos, A. Lozinski, P.Y. Rohan, S. Tomar, Quantifying discretization errors for soft tissue sim- ulation in computer assisted surgery: a preliminary study, Appl. Math. Model. 77 (2020) 709–723.

[22] Z.W. Fang, J.L. Zhang, H.W. Sun, A fast finite volume method for spatial frac- tional diffusion equations on nonuniform meshes, Comput. Math. Appl. 108 (2022) 175–184.

[23] M. Gorakifard, I. Cuesta, C. Saluena, E.K. Far, Acoustic wave propagation and its ap- plication to fluid structure interaction using the cumulant lattice Boltzmann method, Comput. Math. Appl. 87 (2021) 91–106.

[24] M.J. Grote, M. Mehlin, T. Mitkova, Runge-Kutta-based explicit local time-stepping methods for wave propagation, SIAM J. Sci. Comput. 37 (2) (2015) A747–A775.

[25] I. Harari, Dispersion, pollution, and resolution, in: Computational Acoustics of Noise Propagation in Fluids-Finite and Boundary Element Methods, Springer, 2008, pp. 37–56.

[26] D. Hilhorst, M. Vohralík, A posteriori error estimates for combined finite volume- finite element discretizations of reactive transport equations on nonmatching grids, Comput. Methods Appl. Mech. Eng. 200 (5–8) (2011) 597–613.

[27] F. Ihlenburg, I. Babuška, Finite element solution of the Helmholtz equation with high wave number part II: the hp version of the FEM, SIAM J. Numer. Anal. 34 (1) (1997) 315–358.

[28] C. Jansari, J. Videla, S. Natarajan, S.P.A. Bordas, E. Atroshchenko, Adaptive en- riched geometry independent field approximation for 2d time-harmonic acoustics, Comput. Struct. 263 (2022) 106728.

[29] Y. Jin, O.A. González-Estrada, O. Pierard, S.P.A. Bordas, Error-controlled adaptive extended finite element method for 3d linear elastic crack propagation, Comput. Methods Appl. Mech. Eng. 318 (2017) 319–348.

[30] T. Khajah, X. Antoine, S.P.A. Bordas, B-spline fem for time-harmonic acoustic scat- tering and propagation, J. Theor. Comput. Acoust. 27 (03) (2019) 1850059.

[31] T. Khajah, V. Villamizar, Highly accurate acoustic scattering: isogeometric analysis coupled with local high order farfield expansion abc, Comput. Methods Appl. Mech. Eng. 349 (2019) 477–498.

[32] I. Kissami, M. Seaid, F. Benkhaldoun, Numerical assessment of criteria for mesh adaptation in the finite volume solution of shallow water equations, Adv. Appl. Math. Mech. 12 (2) (2020) 503–526.

[33] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, vol. 31, Cambridge University Press, 2002.

[34] R. Luce, B.I. Wohlmuth, A local a posteriori error estimator based on equilibrated fluxes, SIAM J. Numer. Anal. 42 (4) (2004) 1394–1414.

[35] S. Mishra, Ch. Schwab, J. Šukys, Multi-level Monte Carlo finite volume methods for uncertainty quantification of acoustic wave propagation in random heterogeneous layered medium, J. Comput. Phys. 312 (2016) 192–217.

[36] N.M. Newmark, A method of computation for structural dynamics, J. Eng. Mech. Div. 85 (3) (1959) 67–94.

[37] Y. Qu, Z. Li, J. Huang, J. Li, Z. Guan, Elastic full-waveform inversion for surface topography, Geophysics 82 (5) (2017) R269–R285.

[38] S. Schoeder, S. Sticko, G. Kreiss, M. Kronbichler, High-order cut discontinuous Galerkin methods with local time stepping for acoustics, Int. J. Numer. Methods Eng. 121 (13) (2020) 2979–3003.

[39] J. Shragge, B. Tapley, Solving the tensorial 3D acoustic wave equation: a mimetic finite-difference time-domain approach, Geophysics 82 (4) (2017) T183–T196.

[40] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, vol. 22, Springer Science & Business Media, 2013.

[41] E.F. Toro, R.C. Millington, L.A.M. Nejad, Towards Very High Order Godunov Schemes, Springer US, Boston, MA, 2001, pp. 907–940.

[42] R. Velasco-Segura, P.L. Rendón, A finite volume approach for the simulation of nonlinear dissipative acoustic wave propagation, Wave Motion 58 (2015) 180–195.

[43] J. Videla, C. Anitescu, T. Khajah, S.P.A. Bordas, E. Atroshchenko, h- and p-adaptivity driven by recovery and residual-based error estimators for pht-splines applied to time-harmonic acoustics, Comput. Math. Appl. 77 (9) (2019) 2369–2395.

[44] H. Wang, W. Yang, Y. Huang, Adaptive finite element method for the sound wave problems in two kinds of media, Comput. Math. Appl. 79 (3) (2020) 789–801.

[45] X. Wang, B. Hornby, K. Dodds, Dipole sonic response in deviated boreholes pen- etrating an anisotropic formation, in: SEG Technical Program Expanded Abstracts 2002, Society of Exploration Geophysicists, 2002, pp. 360–363.

[46] T. Weinzierl, Higher order methods, in: Principles of Parallel Scientific Computing, Springer, 2021, pp. 249–263.

[47] D. Wu, J. Lv, H. Qian, Parallel domain decomposition schemes based on finite volume element discretization for nonsteady-state diffusion equations on distorted meshes, Comput. Math. Appl. 112 (2022) 97–115.

[48] L.K. Xuan, G.Y. Jin, J.F. Gong, W.P. Zhang, P.J. Ming, Time domain finite volume method for three-dimensional structural-acoustic coupling analysis, Appl. Acoust. 76 (2014) 138–149.

[49] W. Zhang, Y. Zhuang, L. Zhang, A new high-order finite volume method for 3D elas- tic wave simulation on unstructured meshes, J. Comput. Phys. 340 (2017) 534–555.