Improved Approximation of Phase-Space Densities on Triangulated Domains Using Discrete Flow Mapping with p-Refinement

Abstract

We consider the approximation of the phase-space flow of a dynamical system on a triangulated surface using an approach known as Discrete Flow Mapping. Such flows are of interest throughout statistical mechanics, but the focus here is on flows arising from ray tracing approximations of linear wave equations. An orthogonal polynomial basis approximation of the phase-space density is applied in both the position and direction coordinates, in contrast with previous studies where piecewise constant functions have typically been applied for the spatial approximation. In order to improve the tractability of an orthogonal polynomial approximation in both phase-space coordinates, we propose a careful strategy for computing the propagation operator. For the favourable case of a Legendre polynomial basis we show that the integrals in the definition of the propagation operator may be evaluated analytically with respect to position and via a spectrally convergent quadrature rule for the direction coordinate. A generally applicable spectral quadrature scheme for integration with respect to both coordinates is also detailed for completeness. Finally, we provide numerical results that motivate the use of p-refinement in the orthogonal polynomial basis.

1 Introduction

The long-time evolution of dynamical systems is often studied by modelling the propagation of a density distribution $\rho$ in phase-space [16], in particular when the dynamics are ergodic. Conservation of energy (or probability) leads to a phase-space continuity equation known as the generalized Liouville equation

$$\begin{aligned} \frac{\partial \rho }{\partial t} + \nabla \cdot (\mathbf {V}\rho ) = 0, \end{aligned}$$

(1)

where $\mathbf {V}$ is the vector field defining the underlying dynamical system thus

$$\begin{aligned} \dot{X}=\mathbf {V}(X). \end{aligned}$$

(2)

The solutions of this dynamical system $X(\tau )=\phi ^{\tau }(X(0))$ are trajectories in phase-space, where $\phi ^{\tau }$ is the associated flow map. Such models lie at the heart of statistical mechanics and have a wide range of applications including particle tracking in fluids [5, 36], interface tracking in atmospheric sciences [35], the simulation of molecular dynamics [31] as well as illumination or rendering problems in computer graphics [21] or, more generally, the geometrical optics/geometrical acoustics limit of linear wave equations [37].

The method of characteristics leads to an expression for the solution of the generalised Liouville equation in terms of the flow map $\phi ^{\tau }$. A density $\rho$ is transported along trajectories given by $\phi ^{\tau }$ under the action of a linear transfer operator $\mathscr {L}^\tau$ known as the Frobenius–Perron operator [12]

$$\begin{aligned} \rho (X,\tau )=\mathscr {L}^{\tau }[\rho ](X)=\int \delta (X-\phi ^{\tau }(Y))\rho (Y,0)\, \mathrm {d}Y. \end{aligned}$$

(3)

A number of transfer operator based methods have been developed for computing phase-space densities in practical applications. Domain based transfer operator approaches, for example, start by subdividing the phase-space into distinct cells and considering transition rates between these phase-space regions. One of the most common approaches of this type is Ulam’s method [26, 41]. Other methods include wavelet and spectral methods for the infinitesimal Frobenius–Perron operator [15, 20], eigenfunction expansion methods [4] and periodic orbit techniques [12, 27]. For a discussion of convergence properties of the Ulam method see [1, 2]. However, such methods have only found a fairly limited range of applicability, with difficulties arising due to the high-dimensionality of the phase-space.

In this work we consider the case where $\mathbf {V}$ describes a Hamiltonian system, and hence the divergence operator in the generalised Liouville equation (1) reduces to the Poisson bracket as $\nabla \cdot (\mathbf {V}\rho )=\{\rho ,H\}$, where H is the associated Hamiltonian. This class of problems is particularly interesting since it includes, as an important case, short wavelength asymptotic approximations for solving linear wave equations in terms of their underlying ray dynamics. Note that for homogeneous domains with a constant wave speed, and thus a constant and conserved phase-space volume, then the right side of Eq. (3) evaluates to $\rho (\phi ^{-\tau }(X),0)$. In this context, direct “deterministic” methods based on directly tracking swarms of trajectories in phase-space have proved very popular. Such methods are often collectively referred to as ray tracing, see for example [6]. Methods related to ray tracing but tracking the time-dynamics of interfaces in phase-space, such as moment methods and level set methods, have been developed in [14, 28, 32, 43, 44] amongst others. Direct physics-preserving finite difference discretisations of the Liouville equation have also been proposed in Refs. [19, 42]. These methods are best suited to problems with relatively low reflection orders (or problems in one space dimension) and have found applications in acoustics, seismology and optical illumination problems.

Direct ray and interface tracking methods are, however, unsuitable for the modelling of long time dynamics in bounded domains $\varOmega \subset \mathbb {R}^2$. Here, the dynamics (regular or chaotic) depend heavily on the geometry of the boundary of $\varOmega$; multiple reflections of the rays often lead to complicated folding patterns of the associated level-surfaces and an exponential increase in the number of branches that need to be considered. In this regime, the modelling becomes more tractable if one considers the long time stationary solutions to the Liouville equation given by

$$\begin{aligned} \{\rho ,H\}=0. \end{aligned}$$

(4)

In particular, simplified statistical models emerge under additional ergodicity and mixing assumptions on the ray dynamics. A popular tool in the context of vibro-acoustic modelling is Statistical Energy Analysis (SEA) (see for example [29, 30]). A related method for electromagnetic waves is the random coupling model (see [17] and references therein). In SEA, a built-up structure is divided into a set of subsystems and ergodicity of the underlying ray dynamics as well as quasi-equilibrium conditions are imposed. The result is a model where the density in each subsystem is a single degree of freedom, leading to greatly simplified equations based only on flow rates (or coupling loss factors) between subsystems. The assumptions necessary for an SEA treatment are often hard to verify a priori, or are only justified when averaging over an ensemble of similar structures; for example, structures originating from the same production line that are “identical” up to manufacturing tolerances. The limits of SEA have been discussed by Langley and Le Bot amongst others [3, 23,24,25, 38].

In order to extend the statistical approaches described above, we consider numerical schemes for solving the stationary Liouville equation (4) without further assumptions on the nature of the trajectory dynamics. The so-called Dynamical Energy Analysis (DEA) introduced in [37] is based on numerically solving a reformulation of the integral equation (3). In this reformulation, the time-dependent flow map $\phi ^{\tau }$ is replaced with a boundary map $\phi$ that transports trajectories between intersections with the boundary of a finite sub-structure $\varOmega _j\subset \varOmega$, $j=1,\ldots ,N$. Here N can be considered as the number of sub-systems so that

$$\begin{aligned} \varOmega =\cup _{j=1}^N\varOmega _j, \end{aligned}$$

which mirrors an SEA treatment of the problem. Replacing $\phi ^{\tau }$ with $\phi$ simplifies the Frobenius Perron operator (3) to a boundary integral operator, with the resulting advantage of a reduction in the dimensionality. The explicit time-dependence of the density is lost, but the long time stationary behaviour can be recovered by summing over infinitely many iterates of the boundary operator. A trigonometric polynomial basis approximation of the density $\rho$ on the boundary phase-space of each $\varOmega _j$ projects $\rho$ onto finite dimensional space, which corresponds to an SEA model at zeroth order. At higher order the approximation incorporates the geometry dependent ray dynamics, relaxing the underlying ergodicity and quasi-equilibrium assumptions of SEA.

A number of approaches have been suggested in order to obtain efficient discretizations of the boundary operator in DEA. A boundary element method for the stationary Liouville equation (4) was proposed in [9, 10], which extended the DEA approach to larger multi-component structures and three-dimensional applications. It is worth emphasizing here that since the density $\rho$ lives in phase-space, it also needs to be discretized in the direction (or momentum) coordinate. A global orthogonal polynomial basis approximation has typically been used for this purpose. One advantage of a full phase-space formulation such as DEA is that problems due to caustics where ray trajectories focus on a single point in position space are avoided. In phase-space the rays do not intersect since their momentum coordinates are distinct, and it is only after projecting down onto position space that the caustics become apparent; for an example of a DEA simulation including caustics see [11]. Hence caustics do not affect the convergence of the boundary integral model itself, only the post-processing step to compute the density distribution within the solution domain.

A major advantage of DEA is that the choice of subsystem division is no longer critical, since the SEA assumptions related to this criticality have been removed. In fact, one can use the elements of a mesh as substructures in DEA providing huge flexibility and a wealth of potential applications. Applying DEA on the same meshes used for finite element models in industrial applications means that the end-user can access mid-to-high frequencies without the need for extensive remodelling, as would be required for SEA. In Ref. [11] it was shown that local, piecewise constant spatial approximations of the density on the edges of a triangulated surface could be efficiently computed using semi-analytic techniques, which exploit the local geometric simplicity. In this so-called Discrete Flow Mapping (DFM) approach, the trajectory flow is approximated by a discretized flow across a mesh. The technique has since been extended to general convex polygonal mesh elements and to industrial scale applications [8].

Whilst DFM employs orthogonal polynomial approximations of the density in the direction coordinate, the spatial approximation has thus far been limited to piecewise-constant local approximations on the edges of the mesh elements, or subdivisions of those edges, in order to compute the boundary operator using fast semi-analytic methods. This can lead to poor accuracy, particularly in problems with high dissipation where the solutions are dominated by the initial density. In this work, we discuss a p-refinement strategy in both position and direction to gain an improved approximation of the phase-space density on a prescribed mesh. To facilitate the efficient implementation of such a strategy we have derived semi-analytic integration formulae for the projection of the boundary integral operator onto a higher order orthogonal polynomial basis. Here, the integration with respect to the spatial coordinate is performed analytically and the integration with respect to the direction coordinate is done numerically. In fact, a carefully designed strategy for computing the arising multi-dimensional integrals is crucial for the efficient implementation of the method on large meshes. We show in this paper that a direct approach to the analytic spatial integration using the repeated application of the integration by parts rule leads to highly unstable results. However, an indirect approach using the recurrence relations of the Legendre polynomials can be shown to give a stable representation for the integral in terms of modified Bessel functions. Furthermore, an efficient computation of the integrals, independently of the choice of orthogonal polynomial basis functions, is possible using a 2D adaptive and spectrally convergent quadrature method. In this case one can exploit the smoothness of the integrands in the boundary operator over appropriately defined subsets.

The remainder of the paper is structured as follows: in Sect. 2 we detail the boundary version of the Frobenius–Perron operator (3) as a model for propagating an initial phase-space density to give a solution to the stationary Liouville equation (4). In Sect. 3 we describe the finite dimensional approximation of the phase-space density in both the position and direction coordinates. The computational procedure for obtaining the discretized boundary integral operator is described, including three approaches for computing the spatial integrals: an unstable direct formula, a stable recursion formula and a spectral quadrature method. Finally, in Sect. 4 we give a selection of numerical results; we reproduce the well-known decay rate for propagation into free space and perform computations on several bounded domain configurations with varying degrees of symmetry and complexity.

2 Propagating Phase-Space Densities via Integral Operators

In this section we describe the propagation of a phase-space density via integral operators in polygonal multi-domain problems, for example, triangle meshes. We begin by introducing a boundary map $\phi$ and a boundary integral operator $\mathscr {B}$. We then discuss how an initial boundary density can be obtained from a point source and how the final stationary density can be computed from the Neumann series for the operator $\mathscr {B}$.

Fig. 1

Left the boundary map $\phi (s,p_s)$ and the phase-space coordinates $(s,p_s)$ on the boundary $\varGamma$. Right the multi-domain boundary map $\phi _{i,j}(s_j,p_j)$ taking the phase-space coordinates $(s_j,p_j)$ on the boundary $\partial \varOmega _j$ to the phase-space coordinates $(s_i,p_i)$ on $\partial \varOmega _i$

2.1 Boundary Integral Operators

Consider a polygonal domain $\varOmega \subset \mathbb {R}^2$ with boundary $\varGamma$ and the phase-space coordinates $Y_s=(s,p_s)$ on $\varGamma$ as illustrated in the left plot of Fig. 1. Here s parametrizes $\varGamma$ and $p_s$ denotes the direction component tangential to $\varGamma$ at s. Figure 1 also depicts the boundary map $\phi (s,p_s)=(s',p'_{s})$, which takes $(s,p_s)$ to a point $s'\in \varGamma$ with direction $p'_{s}$. We may write this simply as $X_s=\phi (Y_s)$. In general, the form of the boundary map $\phi$ depends on the specific problem. In our considerations, $\phi$ obeys the law of (specular) reflection.

A phase-space density $\rho$ is transported from the phase-space on the boundary to the next boundary intersection via the boundary integral operator [11]

$$\begin{aligned} \mathscr {B}[\rho ](X_s):= \int w(Y_s) \delta (X_s-\phi (Y_s)) \rho (Y_s) \, \mathrm {d}Y_s. \end{aligned}$$

(5)

Note that for $w\equiv 1$, this operator corresponds to the Perron-Frobenius operator (3) for the boundary map $\phi$. The weighting term w is introduced to incorporate absorption factors as well as reflection/transmission coefficients.

In this paper we are concerned with modelling the Hamiltonian ray dynamics of position $\mathbf {r}$ and momentum (slowness vector) $\mathbf {p}$. The Hamiltonian of the trajectory equations is $H=c|\mathbf {p}|=1$, where c is the phase velocity. Conservation of energy defines the integration limits for $p_s\in (-c^{-1},c^{-1})$, since $p_s$ is the tangential component of the direction $\mathbf {p}$. Thus the (double) integral in (5) is over $\varGamma \times (-c^{-1},c^{-1})$. We note that this corresponds to the Hamiltonian for the ray trajectories obtained in the geometrical optics limit for the reduced wave (or Helmholtz) Eq. [14]

$$\begin{aligned} \triangle u + k^2 u =0, \end{aligned}$$

(6)

with $k=\omega /c$ the acoustic wavenumber at angular frequency $\omega$.

The stationary density on the boundary $\varGamma$ induced by a given initial boundary density $\rho _{0}^{\varGamma }$ can be obtained from a Neumann series via

$$\begin{aligned} \rho ^{\varGamma }(s,p_s) = \sum _{n=0}^{\infty }\mathscr {B}^{n}\left[ \rho _{0}^{\varGamma }\right] (s,p_s) = (I-\mathscr {B})^{-1}\left[ \rho _{0}^{\varGamma }\right] (s,p_s), \end{aligned}$$

(7)

where $\mathscr {B}^{n}$ models trajectories undergoing n reflections at the boundary and I is the identity operator. Once the stationary density $\rho ^{\varGamma }$ is found, the interior density $\rho _{\varOmega }$ can then be obtained by projecting down onto position space in $\varOmega$ using [9]

$$\begin{aligned} \rho _{\varOmega }(\mathbf {r}) = \frac{1}{c^2} \int _{\varGamma } \rho ^{\varGamma }(s,p_s) \frac{\cos (\vartheta (s,\mathbf {r}))}{|\mathbf {r}-\mathbf {r}_s|} \, \mathrm {d}s. \end{aligned}$$

(8)

Here $\mathbf {r}\in \varOmega$ is a prescribed solution point and $\mathbf {r}_s$ denotes the Cartesian coordinates of the point $s\in \varGamma$. The angle $\vartheta (s,\mathbf {r})$ is formed between the normal vector to $\varGamma$ at s and the direction vector pointing towards the point $\mathbf {r}$ from the point $s\in \varGamma$.

2.2 Source Terms

In this section we describe how to obtain an initial density $\rho _{0}^{\varGamma }$ from a point source at $\mathbf {r}_0\in \varOmega$ emanating trajectories continuously in all directions. Starting from the high frequency asymptotics for the fundamental solution of the reduced wave Eq. (6) in $\mathbb {R}^2$, $\rho _{0}^{\varGamma }$ is derived from the acoustic energy density

$$\begin{aligned} \varepsilon =\rho ^fk^2|u|^2, \end{aligned}$$

(9)

where u is interpreted as the velocity potential in a fluid medium of density $\rho ^f$. We then obtain a phase-space source density $\rho _0$ by lifting $\varepsilon$ to phase-space with a direction vector $\mathbf {p}_0$, corresponding to the direction coming from $\mathbf {r}_0$. Applying this and replacing u in (9) with the high frequency fundamental solution of (6) yields [9]

$$\begin{aligned} \rho _0(\mathbf {r},\mathbf {p};\mathbf {r}_0) = \rho ^fk^2\left| \frac{i}{4}H_{0}^{1} \left( \frac{\omega |\mathbf {r}-\mathbf {r}_0|}{c} \right) \right| ^2 \delta (\mathbf {p}-\mathbf {p}_0) \underset{\omega \gg 1}{\approx } \rho ^fk\frac{\delta (\mathbf {p}-\mathbf {p}_0)}{8\pi |\mathbf {r}-\mathbf {r}_0|}, \end{aligned}$$

(10)

where $H_{0}^{1}$ is the zeroth order Hankel function of the first kind. In the high frequency limit, the energy density (10) decays in all directions according to $|\mathbf {r}-\mathbf {r}_0|^{-1}$. Considering (10) in the high frequency limit, it can be shown [9] that the initial boundary density distribution from the source point $\mathbf {r}_0$ is

$$\begin{aligned} \rho _{0}^{\varGamma }(s,p_s;\mathbf {r}_0) = \frac{\omega \rho ^f\cos (\vartheta (s,\mathbf {r}_0)) \delta (p_s-p_{s_{0}})}{8\pi |\mathbf {r}_s - \mathbf {r}_0|}, \end{aligned}$$

(11)

where $p_{s_0}$ denotes the component of $\mathbf {p}_0$ tangential to $\varGamma$. The resulting boundary density (11) is equivalent to the source density on the boundary producing the same interior density as the original source distribution (10).

2.3 Multi-domain Problems

A generalization to multi-domain problems such as triangle meshes with sub-domains $\varOmega _{j}$, $j=1,\dots ,N$, is straightforward by introducing a multi-domain boundary map $\phi _{i,j}$ and a weight function $w_{i,j}$ describing the flow from the boundary of the domain $\varOmega _j$ to the boundary of the domain $\varOmega _i$, see the right plot of Fig. 1. Note that $\varOmega =\cup _{j=1}^{N}\varOmega _{j}$ becomes the union of all sub-domains and $\varGamma$ becomes the union of all sub-domain boundaries thus $\varGamma =\cup _{j=1}^{N}\partial \varOmega _{j}$. Each sub-domain $\varOmega _j$ has its own phase-space boundary coordinates $(s_j,p_j)$. We then define the boundary integral operator $\mathscr {B}_{i,j}$, which transports the phase-space density $\rho ^{\varGamma }$ from the boundary phase-space of $\varOmega _j$ to the boundary phase-space of $\varOmega _i$ as illustrated in the right plot of Fig. 1. If the properties of two neighboring domains $\varOmega _j$ and $\varOmega _i$ are different, for example if $c_i\ne c_j$, where $c_i$ is the propagation speed in $\varOmega _i$ (likewise for $c_j$ in $\varOmega _j$), then the weight function $w_{i,j}$ will account for the probability of transmission or reflection at the common edge. Note that the case of transporting the density within a sub-domain is included in this formulation when $i=j$. The operator $\mathscr {B}$ is then constructed from the set of inter-domain operators $\mathscr {B}_{i,j}$. In order to obtain a discrete representation of the phase-space density $\rho ^{\varGamma }$ and the operator $\mathscr {B}_{i,j}$, we consider a finite basis approximation as described in the next section.

3 Discretization

In this section we discuss the discretization of the boundary operator as well as computational issues associated with the efficient and fast implementation of DFM on triangle meshes. We detail how the two-dimensional integral in space and direction can be separated and provide analytical formulae for the spatial integral. We note that a direct approach to evaluating the spatial integral analytically is unstable in general for non-constant spatial basis functions. An alternative stable iterative approach is presented instead and compared with an efficient, and widely applicable, spectrally convergent quadrature strategy. We conclude this section by outlining a computational algorithm for DFM with p-refinement on triangle meshes.

The finite dimensional representation of $\rho ^{\varGamma }$ is given by a basis approximation of the form

$$\begin{aligned} \rho ^{\varGamma }(s_j,p_j)\approx \sum _{l=1}^{N_j}\sum _{m=0}^{N_s}\sum _{n=0}^{N_p}\rho _{(j,l,m,n)} \hat{P}_{m}^{l}(s_j) \tilde{P}_{n}(p_j), \qquad j=1,\dots ,N, \end{aligned}$$

(12)

where $N_j$ is the number of the boundary elements for $\varOmega _j$, $N_{s}$ is the order of the basis expansion in space and $N_{p}$ is the order of the basis expansion in direction. Note that the values of $N_s$ and $N_p$ may vary depending on the sub-domain $\varOmega _j$. However, for simplicity in the sequel we only consider triangle meshes with one boundary element per edge ($N_j\equiv {3}$ for all j) and we take constant values for $N_s$ and $N_p$. Refinement of the spatial approximation can therefore take place via both increasing $N_s$ and refining the triangle mesh. The summand in (12) comprises a product of the unknown expansion coefficients $\rho _{(j,l,m,n)}$, the basis functions $\hat{P}_m^{l}$ in the spatial coordinate and the basis functions $\tilde{P}_{n}$ in the direction coordinate. We assume that both sets of basis functions (space and direction) are orthonormal, and that $\hat{P}_m^{l}$ is smooth within the $l{\text {th}}$ boundary element and $\tilde{P}_{n}$ is smooth on $(-c^{-1},c^{-1})$.

In this work we concentrate on the case when $\hat{P}_m^{l}$ is given by a scaled Legendre polynomial of the form

$$\begin{aligned} \hat{P}_m^{l}(s_j) = \sqrt{\frac{2}{A_{l,j}}} P_m(t)\chi _{l}(s_j), \quad t=\frac{2s_j-A_{l,j}-2\sum _{\kappa <l}A_{\kappa ,j}}{A_{l,j}} \in [-1,1), \end{aligned}$$

where $t\in [-1,1)$ is a parametrization of the $l{\text {th}}$ boundary element of $\varOmega _j$, with length $A_{l,j}$, defined here as a function of $s_j\in [\sum _{\kappa<l}A_{\kappa ,j},A_{l,j}+\sum _{\kappa <l}A_{\kappa ,j})$. Also, $P_m$ denotes the Legendre polynomial of degree m and $\chi _{l}(s_j)$ is an indicator function taking the value 1 when $s_j$ is in the $l{\text {th}}$ element of $\varOmega _j$, and zero otherwise. We likewise focus on the case when $\tilde{P}_{n}$ is given by the scaled Legendre polynomial:

$$\begin{aligned} \tilde{P}_{n}(p_j)= \sqrt{c_j} P_{n}(c_j p_j). \end{aligned}$$

(13)

Note that the Legendre polynomials represent a particularly convenient choice due to their simple orthogonality relation, their spectral convergence for sufficiently smooth approximants and their ease of calculation using a recurrence relation. Later we shall see that this choice also permits highly efficient semi-analytic integration routine. However, we will also describe a fully numerical and spectrally convergent quadrature algorithm for any smooth orthonormal basis choice in order to demonstrate the broader scope and applicability of the proposed methodology.

A Galerkin projection of the operator $\mathscr {B}_{i,j}$ on to the orthonormal basis in the expansion (12) leads to a matrix representation $\mathbf {B}$ with entries given by

$$\begin{aligned} \mathbf {B}_{I,J}=\alpha _{m',n'} \int _{Q_j} w_{i,j}(Y_{j}) \tilde{P}_{n'}(\phi _{p}(Y_{j})) \hat{P}_{m'}^{l'}(\phi _{s}(Y_{j})) \tilde{P}_{n}(p_{j}) \hat{P}_{m}^{l}(s_j) \,\mathrm {d}Y_j, \end{aligned}$$

(14)

where I and J denote the multi-indices $I=(i,l',m',n')$ and $J=(j,l,m,n)$, respectively. The factor $\alpha _{m',n'}$ is a scaling required for orthonormality; for a Legendre polynomial basis then $\alpha _{m',n'}=(2m'+1)(2n'+1)/4$. We have introduced the notation $Y_j=(s_j,p_j)\in Q_j$, where $Q_j=\partial \varOmega _j\times (-c_j^{-1},c_{j}^{-1})$, and $\phi =(\phi _s,\phi _p)$ has been used to separate the boundary map into spatial and directional components. Since the expansion (12) represents a local density approximation and the transfer of energy is only between connected sub-domains, $\mathbf {B}$ is a block-sparse matrix for large N. The coefficients of the expansion (12) can be found by solving the linear system ${\varvec{\rho }} = (\mathbf {I}-\mathbf {B})^{-1}{\varvec{\rho }}_0$, which corresponds to the discretized form of Eq. (7). Here $\varvec{\rho }_0$ and $\varvec{\rho }$ represent the coefficients of the expansions of $\rho _{0}^{\varGamma }$ and $\rho ^{\varGamma }$, respectively, when projected onto the finite dimensional basis. Note that for $\mathbf {I}-\mathbf {B}$ to be invertible, absorption/dissipation must be included in the weights $w_{i,j}$. Once $\varvec{\rho }$ has been computed and substituted into (12), then the interior density $\rho _\varOmega$ can be computed using (8).

3.1 DFM with p-Refinement on Triangle Meshes

We now describe a strategy for computing the entries of the transfer matrix $\mathbf {B}$, given in (14), on the elements of a triangle mesh. Note that the integration over the direction coordinate $p_j$ can instead be performed with respect to the angle $\theta _j$ between the normal vector to $\varGamma$ at $s_j$ and the outgoing trajectory direction using the relation $p_j=c_j^{-1}\sin (\theta _j)$. We will assume that the weight function $w_{i,j}$ can be written in the form

$$\begin{aligned} w_{i,j}(Y_j)=e^{-\mu d(s'_i,s_j)}\lambda _{i,j}(\theta _j), \end{aligned}$$

where the factor $\mu$ describes the rate of dissipation along the direction of propagation. Here we write $\phi _{s}(Y_{j})=s_i'$ for brevity and $d(s'_i,s_j)$ is the Euclidean distance between the points represented by $s_i'$ and $s_j$. The function $\lambda _{i,j}$ gives the reflection/transmission probability, which may be derived from the underlying wave equation by assuming continuity of the pressure and velocity for an incident plane wave [13]. Writing the specific acoustic impedance in $\varOmega _i$ as $z_i=\rho ^f_{i}c_i$, the product of the fluid density and the propagation speed inside $\varOmega _i$, then one obtains the transmission probability

$$\begin{aligned} \lambda _{t}(\theta _j) = \frac{4z_i z_j\cos (\theta _r)\cos (\theta _t)}{(z_i\cos (\theta _r)+z_j\cos (\theta _t))^2}. \end{aligned}$$

(15)

The reflection/transmission probability function is then defined as $\lambda _{i,j}(\theta _j)=\lambda _{t}(\theta _j)$ for $i\ne j$, and $\lambda _{j,j}(\theta _j)=1-\lambda _{t}(\theta _j)$. Note that in (15) above, $\theta _r$ and $\theta _t$ are the $\theta _j$–dependent reflection and transmission angles, respectively; see the left plot of Fig. 2.

In what follows we denote $s_j$ as the coordinate on the initial edge, we label quantities associated with the edge connected to the initial edge in the counter-clockwise direction by ‘$+$’ and we label quantities associated to the third edge by ‘−’. We find

$$\begin{aligned} \theta _{r}^{\pm }(\theta _j) = \pm \varphi _{\pm } - \theta _j, \end{aligned}$$

where $\varphi _{\pm }$ are the interior angles of $\varOmega _j$ between the edge containing $s_j$ and the two possible destination edges, as shown in the left plot of Fig. 2. The transmission angles $\theta _t^{\pm }$ are obtained from the reflection angles $\theta _{r}^{\pm }$ using Snell’s law

$$\begin{aligned} \sin (\theta _{t}^{\pm })= - \frac{c_{i}}{c_{j}}\sin (\theta _{r}^{\pm }), \end{aligned}$$

(16)

where the minus sign comes from the fact that $\theta _t^{\pm }\in (-\pi /2,\pi /2)$ is defined in the local coordinates of the neighboring triangle $\varOmega _i$.

Fig. 2

Left the evolution of trajectories through a triangulated mesh, including reflection and transmission angles $\theta _r$ and $\theta _t$, respectively. Right admissible ranges for the spatial integration defined using $s_j^*$

For $c_i>c_j$, the transmission angle $\theta _t^{\pm }$ in Snell’s law (16) is only defined for $\theta _j$ in the interior of the set $\varTheta$ where

$$\begin{aligned} \varTheta :=[-\pi /2,\pi /2]\cap (\pm \varphi _{\pm }-\arcsin (c_j/c_i),\pm \varphi _{\pm }+\arcsin (c_j/c_i)). \end{aligned}$$

Total internal reflection ($\lambda _t=0$) takes place for $\theta _j\notin \varTheta$, and this should be taken into account when integrating $\lambda _t$ over $\theta _j$. The convergence of numerical quadrature rules can be adversely affected by the branch points of $\arcsin$ at $\pm 1$ in the Snell’s law expression for $\theta _t$ (16) when $\varTheta ^c:=[-\pi /2,\pi /2]\setminus \varTheta$ is non-empty. Furthermore, separating the integral with respect to $\theta _j$ into a sum of integrals over $\varTheta$ and $\varTheta ^c$ is necessary, but not sufficient to address this issue. For integration over $\varTheta ^c$ we may simply set $\lambda _t=0$, but since the branch points are located at the boundary between $\varTheta$ and $\varTheta ^c$, problems still arise for the integration over $\varTheta$. Note that from (16) we obtain

$$\begin{aligned} \cos (\theta _{t}^{\pm })=\sqrt{1-(c_i/c_j)^2\sin ^2(\theta _r^{\pm })} \end{aligned}$$

(17)

and hence the branch points in formula (15) are located at zeros of the square root in (17). A change of variables from $\theta _j$ to the transmission angle $\theta _t^{\pm }$ can be employed to avoid these branch points. Explicitly, we use that

$$\begin{aligned} \theta _j&= \pm \varphi _{\pm } - \theta _r^{\pm }\\&= \pm \varphi _{\pm } + \arcsin ((c_j/c_i)\sin (\theta _{t}^{\pm })), \end{aligned}$$

and $\cos ^2(\theta _r^{\pm })={1-(c_j/c_i)^2\sin ^2(\theta _t^{\pm })}$ to derive

$$\begin{aligned} \mathrm {d}\theta _j&=\frac{c_j\cos (\theta _t^{\pm })}{c_i\cos (\theta _r^{\pm })} \mathrm {d}\theta _t^{\pm }. \end{aligned}$$

Thus

$$\begin{aligned} \lambda _t(\theta _j) \mathrm {d}\theta _j = \frac{4z_iz_j(c_j/c_i)\cos ^2(\theta _t^{\pm })}{\left( z_i\sqrt{1-(c_j/c_i)^2\sin ^2(\theta _t^{\pm })} + z_j\cos (\theta _t^{\pm }) \right) ^2} \mathrm {d} \theta _t^{\pm } =: \lambda (\theta _t^{\pm }) \mathrm {d}\theta _t^{\pm }, \end{aligned}$$

where the term under the square root in the denominator of $\lambda (\theta _t^{\pm })$ is strictly positive for any value of $\theta _t^{\pm }$ since $c_i>c_j$. The integration limits for the new integration variable $\theta _t^{\pm }$ can then be found by considering the minimal and maximal $\theta _j\in \varTheta$.

We now consider the integration with respect to position $s_j$, which will be performed before integrating over $\theta _j$ and hence we integrate over a direction dependent domain of the form $(s_{min}(\theta _j),s_{max}(\theta _j))$. Consider edge number $l\in \{1,2,3\}$ of the triangle $\varOmega _j$ with length $A_{l,j}$, represented by the edge AB in the right plot of Fig. 2. Then $s_j\in [\sum _{\kappa<l}A_{\kappa ,j},A_{l,j}+\sum _{\kappa <l}A_{\kappa ,j}):=[s_A,s_B)$ on that edge and the admissible ranges for the spatial integration are given by

$$\begin{aligned} s_{min}^{-}(\theta _j) = s_A, \quad s_{max}^{-}(\theta _j) = s^*_j(\theta _j) \end{aligned}$$

(18)

for trajectories transported to the − edge, and

$$\begin{aligned} s_{min}^{+}(\theta _j) = s^*_j(\theta _j), \quad s_{max}^{+}(\theta _j) = s_B, \end{aligned}$$

(19)

for trajectories transported to the $+$ edge. We define the position $s^*_j(\theta _j)$ to be the starting point of the trajectory on the $l{\text {th}}$ edge, which is mapped to the vertex of $\varOmega _j$ opposite to the edge l (shown as vertex C in Fig. 2). We also enforce that

$$\begin{aligned} \mathop {\min }\limits _{\theta _j}(s_j^*)=s_A\,\,\,\ \mathrm {and}\,\,\, \mathop {\max }\limits _{\theta _j}(s_j^*)=s_B. \end{aligned}$$

The behaviour of $s_j^*$ as a function of $\theta _j$ will be important when considering the integral over $\theta _j$. In particular, note that for $\theta _j\leqslant -(\pi /2)+\varphi _+$ then $s_j^*=s_B$ and trajectories only map to the − edge. Likewise, for $\theta _j\geqslant (\pi /2)-\varphi _-$, then $s_j^*=s_A$ and trajectories only map to the $+$ edge. For the intermediate case, $-(\pi /2)+\varphi _+<\theta _j<(\pi /2)-\varphi _-$, then $s_j^*$ is a linear function of $\theta _j$ and takes values along the interior of the edge l. Owing to these abrupt changes from constant to linear behavior, $s_j^*(\theta _j)$ is not smooth at the points $\theta _j=-(\pi /2)+\varphi _+$ and $\theta _j=(\pi /2)-\varphi _-$. This leads to a natural subdivision of the range of integration for $\theta _j$ given by

$$\begin{aligned} \left( -\frac{\pi }{2},\frac{\pi }{2}\right) =\left( -\frac{\pi }{2},-\frac{\pi }{2}+ \varphi _+\right] \cup \left( -\frac{\pi }{2}+\varphi _+,\frac{\pi }{2}-\varphi _ -\right) \cup \left[ \frac{\pi }{2}-\varphi _-,\frac{\pi }{2}\right) \end{aligned}$$

(20)

and a corresponding splitting for the numerical integration over $\theta _j$. Note that the integral over either the first or last of these sub-intervals (20) will always be equal to zero, depending on whether the propagation is to the $+$ or − edge, respectively, since the corresponding spatial integration will have equal lower and upper limits. Note that these subdivisions are always applied, in contrast with the subdivisions at points of total internal reflection discussed before, which are only applied when $c_i>c_j$.

In summary, the entries of the matrix $\mathbf {B}$ (14) can be rewritten in the form:

$$\begin{aligned} \begin{aligned} \mathbf {B}_{I,J} =&\frac{\alpha _{m',n'}}{c_j}\int _{-\pi /2}^{\pi /2} \lambda _{i,j}(\theta _j) \tilde{P}_{n'}(p'_{i}) \tilde{P}_{n}(p_{j}) \cos (\theta _j)\\&\times \left[ \int _{s_{min}(\theta _{j})}^{s_{max}(\theta _{j})} \hat{P}_{m'}^{l'}(s'_{i}) \hat{P}_{m}^{l}(s_{j}) e^{-\mu d(s'_i,s_j)} \, \mathrm {d}s_{j}\right] \, \mathrm {d}\theta _{j}, \end{aligned} \end{aligned}$$

(21)

where $p_i'=c_i^{-1}\sin (\theta _i')$ and $\theta _i'=\theta _r$ if $i=j$, or $\theta _i'=\theta _t$ otherwise. The distance $d(s'_i,s_j)$ and the value of $s'_i$ both depend on whether the ray hits the previous or the next neighboring edge (oriented counter-clockwise) independently whether $i=j$ or not. In fact, in any convex polygonal sub-domain $\varOmega _j$, one can show via geometrical arguments that both $d(s'_i,s_j)$ and $s'_i$ are linear functions of $s_j$. After all subdivisions we re-scale each of the sub-integrals over $\theta _j$ and $s_j$, required for computing (21), and write them in the following general form

$$\begin{aligned} \mathfrak {I} \propto \int _{-1}^{1} f(\theta ) \left[ \int _{-1}^{1} g_{m,m'}(s,\theta ) e^{a(\theta )s+b(\theta )} \, \mathrm {d}s \right] \, \mathrm {d}\theta = \int _{-1}^{1} f(\theta ) \mathfrak {I}_{m,m'}(\theta )\, \mathrm {d}\theta . \end{aligned}$$

(22)

In principle, we could now compute $\mathfrak {I}$ using any quadrature method, and since care has been taken to preserve smoothness of the integrand, then a spectrally convergent method would be preferable. The efficiency of computing $\mathfrak {I}$ in (22) can be further increased if $\mathfrak {I}_{m,m'}(\theta )$ can be evaluated analytically as suggested for the case of piecewise constant spatial basis functions in [11]. We discuss this in the following section for the case when the spatial basis functions are chosen as scaled Legendre polynomials.

3.2 Exact Spatial Integration

In this section we present two formulae that can be used to evaluate the spatial integral $\mathfrak {I}_{m,m'}(\theta )$ in (22) exactly when $\hat{P}_{m}^{l}(s_j)$ are given by scaled Legendre polynomials and we have prescribed a fixed value of $\theta$. Hence $a(\theta )$ and $b(\theta )$ may be considered as constants and $g_{m,m'}(s,\theta )$ is a product of Legendre polynomials of the form $g_{m,m'}(s,\theta )=P_{m}(a_1s+b_1)P_{m'}(a_2s+b_2)$. Henceforth, we omit $\theta$ from the notation when referring to these quantities, writing them as a, b and g(s). The integral $\mathfrak {I}_{m,m'}(\theta )$ may be computed directly using integration by parts as

$$\begin{aligned} \begin{aligned} \mathfrak {I}_{m,m'}(\theta )&= \int _{-1}^{1} g(s) e^{as+b} \, \mathrm {d}s = \left. \left[ \sum _{i=0}^{m+l} \frac{(-1)^i}{a^{i+1}} g^{(i)}(s)\right] e^{as+b} \right| _{-1}^{1} \\&= \left. \left[ \sum _{i=0}^{m+l} \frac{(-1)^i}{a^{i+1}} \sum _{j=0}^{i} \begin{pmatrix} i\\ j \end{pmatrix} a_1^{i-j} a_2^{j} P_{m}^{(i-j)}(a_1s+b_1) P_{m'}^{(j)}(a_2s+b_2) \right] e^{as+b} \right| _{-1}^{1}, \end{aligned} \end{aligned}$$

(23)

where the superscript (i) denotes the $i{\text {th}}$ derivative. Note that the spatial derivatives of the Legendre polynomials can be evaluated via a recursive relationship, and that for $m=m'=0$, the integral has the relatively simple form $\mathfrak {I}_{0,0}(\theta ) =(e^{as+b}/a)|_{-1}^{1}$ as used in [11]. Unfortunately, the formula (23) is numerically unstable for general m and $m'$ (see Fig. 3), since the coefficients in the expansion grow with each derivative. Thus the total sum is highly sensitive to round-off errors. In addition, the coefficients a, b, $a_1$, $b_1$, $a_2$ and $b_2$ are geometry and dissipation rate dependent meaning that (23) may be stable in some mesh elements, but unstable in others.

Fig. 3

Comparison between direct, recursive and numerical methods for evaluating $\mathfrak {I}_{m,m'}(\theta )$ with $\theta \in (-1,1)$. Left integral values for polynomial degrees $m=2$ and $m'=0$. Right integral values for polynomial degrees $m=m'=2$. Parameter values: $c=1$ and $\mu =0.005$

An alternative approach is to derive a recursive formula using the properties of $g(s)=P_{m}(a_1s+b_1)P_{m'}(a_2s+b_2)$. In particular, one can show for each fixed value of $\theta$ that either $|a_1|=1$ and $b_1=0$, or $|a_2|=1$ and $b_2=0$. Assuming without loss of generality that $b_1=0$ and re-scaling via $\bar{s}=a_1s$ we find

$$\begin{aligned} \mathfrak {I}_{m,m'}(\theta ) = \int _{-1}^{1} e^{as+b} P_{m}(a_1 s) P_{m'}(a_2s+b_2) \, \mathrm {d} s = \frac{1}{a_1} \int _{-1}^{1} e^{\tilde{s}} P_{m}(\bar{s}) P_{m'}(\hat{s}) \, \mathrm {d} \bar{s}, \end{aligned}$$

where $\tilde{s}=(a/a_1)\bar{s}+b$ and $\hat{s}=(a_2/a_1)\bar{s}+b_2$. The well-known Legendre polynomial recurrence relationship for $P_{m'}(\hat{s})$ is given by

$$\begin{aligned} (m'+1)P_{m'+1}(\hat{s})=(2m'+1)\hat{s}P_{m'}(\hat{s})-m'P_{m'-1}(\hat{s}). \end{aligned}$$

(24)

Applying this relation together with the definition of $\hat{s}$ in the right-hand expression for $\mathfrak {I}_{m,m'}$ above yields

$$\begin{aligned} \mathfrak {I}_{m,m'} = \frac{a_2(2m'-1)}{a_1^2 m'} \int _{-1}^{1} e^{\tilde{s}} \bar{s} P_{m}(\bar{s}) P_{m'-1}(\hat{s}) \, \mathrm {d} \bar{s} + \frac{b_2(2 m'-1)}{m'} \mathfrak {I}_{m, m'-1} - \frac{m'-1}{m'} \mathfrak {I}_{m,m'-2}. \end{aligned}$$

Applying the Legendre polynomial recurrence relationship once more to $\bar{s}P_{m}(\bar{s})$, we obtain a recurrence relationship for $\mathfrak {I}_{m,m'}(\theta )$ with $m,m'=0,1,2,\dots$ as

$$\begin{aligned} \begin{aligned} \mathfrak {I}_{m,m'} =&\, \frac{a_2(m+1)(2m'-1)}{a_1 (2m+1) m'}\mathfrak {I}_{m+1,m'-1} + \frac{a_2 m (2m'-1)}{a_1(2m+1)m'}\mathfrak {I}_{m-1,m'-1} \\&\quad + \frac{b_2(2m'-1)}{m'} \mathfrak {I}_{m,m'-1} - \frac{m'-1}{m'} \mathfrak {I}_{m,m'-2}. \end{aligned} \end{aligned}$$

(25)

The above recurrence formula may be initiated using the following integral relation:

$$\begin{aligned} \mathfrak {I}_{m,0} =\frac{e^b}{a_1}\int _{-1}^{1} e^{(a/a_1)\bar{s}} P_{m}(\bar{s}) \, \mathrm {d} \bar{s} = \frac{e^b}{a_1}\sqrt{\frac{2\pi a_1}{a}} I_{m+1/2}\left( \frac{a}{a_1}\right) , \end{aligned}$$

which can be derived from ([18], Eq. 7.243–5) and where $I_{m+1/2}$ is the modified Bessel function of the first kind. The computational implementation in the present work has been performed in MATLAB, where we make use of the existing functions for computing the modified Bessel functions. Alternatively, the values of $\mathfrak {I}_{m,0}$ can be computed using the recursion relation for the Legendre polynomial derivatives. Unfortunately, the resulting recursion relation for $\mathfrak {I}_{m,0}$ is unstable in the forward direction. However, it may be performed in the backward direction using continued fractions as detailed in Sect. 6.6.2 of Ref. [34].

We illustrate the numerical instability of the direct formula (23) by computing $\mathfrak {I}_{m,m'}(\theta )$ for $\theta \in (-1,1)$ on a unit sided equilateral triangle with a specified initial edge (due to symmetry the result will be independent of this choice). In Fig. 3, the direct formula (23) is compared against the recursive solution (25) and a numerical solution using Clenshaw–Curtis quadrature [39] for different combinations of the polynomial degrees m and $m'$. In both examples we set the dissipation rate to be $\mu =0.005$ and the propagation speed as $c=1$.

The left plot of Fig. 3 shows the case $m=2$ and $m'=0$, where all three approaches compute the same values for the integral to double precision. The numerical integration was performed using three integration subregions $\theta _j\in (-\pi /2,-\pi /6]\cup (-\pi /6,\pi /6)\cup [\pi /6,\pi /2)$ as described in (20) since $\varphi _{\pm }=\pi /3$. The lack of smoothness of the integral $\mathfrak {I}_{m,m'}(\theta )$ at the numerical integration subdivision points $\theta =2\theta _j/\pi =\pm 1/3$ is also clear in Fig. 3, emphasizing the motivation for performing this subdivision. The right plot of Fig. 3 shows the case when $m=m'=2$, which gives rise to large deviations between the direct formula and other two methods. Note that, in general, the discrepancy increases when the degree of the basis functions is increased and/or the dissipation factor $\mu$ is decreased. The right plot of Fig. 3 shows that the direct formula (23) is already numerically unstable for relatively low degree basis functions. For the MATLAB implementation here, we find that using the recursive formula (25) is typically between two and three times faster than applying the spectral quadrature method. In order to demonstrate the flexibility of the method in accommodating general smooth orthonormal basis approximations in (12), we propose an efficient adaptive quadrature strategy in the following section.

3.3 Efficient Adaptive Quadrature Rule

In this section we describe the spectrally converging adaptive quadrature strategy that we use to compute the integral entries of the matrix $\mathbf {B}$ (21) for general smooth orthonormal basis functions $\hat{P}_{m}^{l}(s_j)$ and $\tilde{P}_{n}(p_j)$. Note that whenever possible, it is preferable to adopt analytic spatial integration strategies for reasons of computational efficiency. However, in order to demonstrate the wider applicability of the proposed methodology, and for completeness, we detail a rapidly converging fully numerical strategy. We propose an adaptive Clenshaw–Curtis method [39, 40] and note that spectral convergence property depends crucially on the regularity of the integrand. Preservation of the rapid convergence property is the main motivation for the integral subdivision strategy discussed in Sect. 3.1.

The subdivision of the integral with respect to $\theta _j$ ensures that integrands are smooth, or even analytic, in a large region of the complex plane. This leads to spectral convergence with a low number of quadrature points in both dimensions, which means that pre-computing and storing the Clenshaw–Curtis quadrature nodes and weights is feasible. In addition, Clenshaw–Curtis quadrature is progressive since the set of nodes for a $(2M+1)$-point quadrature rule includes all the nodes of an $(M+1)$-point rule allowing the repeated use of previously computed data. The result from a $(2M+1)$-point rule is compared with the result from an $(M+1)$-point rule and the code terminates when a desirable tolerance level is achieved. The complete algorithm for the computation of a single entry of the matrix $\mathbf {B}$ (21) with fixed $I=(i,l',m',n')$ and $J=(j,l,m,n)$ is as follows.

• 1. Let $\alpha =0$ and fix $M_0+1$ as an initial number of scaled Clenshaw–Curtis quadrature nodes. Define a prescribed tolerance level "tol" for the accuracy of any quadrature rules.

• 2. Calculate the $M_\alpha +1$ values of $\theta _j$ prescribed by the Clenshaw–Curtis quadrature nodes, then perform the following operations for each value of $\theta _j$:

  • compute the reflection/transmission angle $\theta '_i$ as described in the text immediately after Eq. (21);

  • compute the basis functions $\tilde{P}_{n}(p_j)$ and $\tilde{P}_{n'}(p'_i)$. For a Legendre polynomial basis this can be done using Eq. (13) and the recurrence (24);

  • compute the probability function $\lambda _{i,j}(\theta _j)$ as described in the text immediately after Eq. (15);

  • find the admissible ranges $(s_{min}(\theta _j),s_{max}(\theta _j))$ for the integration with respect to $s_j$ using Eqs. (18) and (19).

• 3. If $\hat{P}_{m}^{l}(s_j)$ are scaled Legendre polynomials, then compute the spatial integral analytically as described in Sect. 3.2. Otherwise:

  • (a) Let $\beta =0$ and fix $m_0+1$ as an initial number of scaled Clenshaw–Curtis quadrature nodes.

  • (b) Calculate the $m_\beta +1$ values of $s_j$ prescribed by the Clenshaw–Curtis quadrature nodes, then perform the following operations for each value of $s_j$:

    • compute the point $s'_i$ and the distance $d(s_i',s_j)$ using the geometry of the mesh and the vector prescribed by $(s_j,p_j)$. On a triangle mesh one can directly compute d as a linear function of $s_j$ using the sine rule;

    • compute the orthogonal polynomial basis functions $\hat{P}_{m}^{l}(s_j)$ and $\hat{P}_{m'}^{l'}(s'_{i})$. The computation will depend on the chosen basis.

  • (c) Compute the integral over $s_j$ appearing in Eq. (21) using $m_\beta$ quadrature nodes. If the tolerance "tol" has not been reached, then add one to $\beta$, let $m_{\beta }=2^\beta m_0$ and return to step (b). If "tol" has been reached then proceed to step 4.

• 4. Compute the integral over $\theta _j$ appearing in Eq. (21) using $M_\alpha$ quadrature nodes. If the tolerance "tol" has not been reached, then add one to $\alpha$, let $M_{\alpha }=2^\alpha M_0$ and return to step 2. If "tol" has been reached then use the computed value as an entry of the matrix $\mathbf {B}$.

Note that for a Legendre polynomial basis, the best strategy is to pre-compute the basis functions up to the maximum order considered whenever possible due to the use of a recurrence formula. The relevant steps in the algorithm above would then be replaced by a call to the pre-computed values whenever these are available. The spectral convergence of the proposed quadrature scheme is demonstrated in “Appendix”. We apply the above steps with "tol"$=10^{-12}$ to obtain the numerical results in the next section.

4 Numerical Results

In this section we primarily study numerical examples which demonstrate the improvement in the phase-space density approximation due to p-refinement in space, since this is the new contribution of this work. We have studied the effects of mesh refinement and higher order approximations in the direction basis in our initial work on DFM [11]. We first study the convergence of DFM for a free space propagation problem, which motivates the use of higher order basis approximations in space. We then consider the problem of computing the energy density in the vicinity of a point source for a dissipative system, including the effect of higher order spatial approximations and local mesh refinement. A closed cavity with a point source at different locations is then considered, and as a final example we model a coupled two-cavity system with piecewise constant propagation speeds. All meshes shown in the results that follow were generated using the DistMesh mesh generator for MATLAB [33].

Fig. 4

Left energy propagation from the point source at the origin with $N_s=3$ and $N_p=300$. Right the absolute error of the energy density averaged over annular segments of width $\varDelta _R=1$ (indicated by the white dashed lines in the left plot). Parameter values for both sub-plots: $\omega =100\pi$, $c=1$, $\rho ^f=1$ and $\mu =0$

4.1 Free Space Propagation

The ray tracing approximation of the energy density $\rho$ emanating from a point source at $\mathbf {r}=\mathbf {r}_0$ in two-dimensional free-space may be obtained from (10) as

$$\begin{aligned} \rho (R)=\frac{\rho ^fk}{8\pi R}, \end{aligned}$$

(26)

where $R=|\mathbf {r}-\mathbf {r}_0|$. Discrete flow mapping is a smoothed ray tracing algorithm in the same way as DEA and, in general, only retrieves the full ray-tracing solution in the limit where $N_s$ and $N_p$ tend to infinity [37]. Figure 4 indicates the convergence of DFM to the ray tracing solution (26), where the triangle mesh has complete transmission ($\lambda _t=1$) across all edges in order to replicate free-space. In these computations we have chosen parameters $\omega =100\pi$, $c=1$, $\rho ^f=1$ and $\mu =0$. The plotted values are computed at the centroid of each triangle and are presented on a logarithmic scale.

In the right plot of Fig. 4 we show the averaged relative error for the energy density computed by DFM with different orders of approximation in space, and a fixed order in direction. The relative error values are averaged over annular segments with radial width $\varDelta _R=1$, as indicated by the white dashed lines in the left plot of Fig. 4. Note that we have excluded the average error in the circular segment containing the source point since the DFM result here is just the direct (and exact) source contribution. Clearly, as the order of the spatial approximation $N_s$ increases, the error is correspondingly decreasing. The approximation in direction is fixed with a relatively large choice of $N_p=300$, since the free-space problem only gives propagation in a single direction, radially outward from the source. As such, the analytic solution in phase-space will behave as a delta-distribution in the direction coordinate, making it challenging to model using a polynomial basis. An identical value for the direction basis order was used in all cases for consistency, and to clearly show the convergence related to increasing the spatial approximation order only.

The estimated order of convergence (EOC) when the number of degrees of freedom is doubled by changing from $N_s=0$ to $N_s=1$ ranges between 1.8 in the closest annular segment to the source point and around 3.5 in the furthest away two segments. Here the EOC is computed via

$$\begin{aligned} \text {EOC}(N_s)=\log _{(N_s+2)/(N_s+1)}\left( \frac{\text {Error}(N_s+1)}{\text {Error}(N_s)}\right) \end{aligned}$$

for $N_s=0$. For the increase from 2 to 3 degrees of freedom per element edge ($N_s=1$ to $N_s=2$), we compute EOC(1), which increases compared with EOC(0) to around 3.5 close to the source and to 5 in the most far field sector. Finally, for the increase from $N_s=2$ to $N_s=3$, EOC(2) does not simply increase further from the source as before. In the closest sector to the source the EOC has reduced to 1.5 and likewise, in the two furthest away sectors we see a similar reduction in the EOC to around 3. However, in the intermediate sectors between $R=2$ and $R=4$ we see an increase in the EOC to as high as 8.4 for the segment from $R=2$ to $R=3$. The saturation of the error at around 3e−4 when $N_s=3$ suggests that here we are close to the maximum accuracy possible on this triangulation with $N_p=300$.

The example shown in Fig. 4 has a relatively small and simple mesh. A basis order of $N_p=300$ would be computationally infeasible for complex built-up engineering applications. However, the directional dependence of the energy density in complex structures is often relatively smooth due to multiple reflections and irregular geometry, which makes such a large choice of $N_p$ unnecessary. In general, engineering applications will also require the study of dissipative problems and hence we now change the dissipation factor to be $\mu =1$ in the free-space radiation problem considered above. Figure 5 shows the numerical solution for this problem with $N_p=12$ and with either a constant or a quadratic spatial basis.

The lower right plot of Fig. 5 shows the ray tracing solution (26) with an additional dissipation factor for consistency with the DFM calculation. The energy density predicted by a DFM simulation with $N_s=0$ and $N_p=12$ is shown in the upper left plot. Spurious localization and shadowing effects arise due to the piecewise constant approximation in space introducing discontinuities at the vertices of the mesh. The upper right plot of Fig. 5 shows another DFM simulation on the same mesh, but here with $N_s=2$ and $N_p=12$. The circular symmetry of the energy density shown in the exact solution has been preserved showing a distinct qualitative improvement from the piecewise constant approximation. The higher order simulation leads to an approximate doubling of the computation time, but even so, the simulations presented here took less than 30 s.

The shadowing effects produced in the upper left plot of Fig. 5 can also be reduced by applying local mesh refinement near the source point. The lower left plot of Fig. 5 shows that the local mesh refinement strategy, which was implemented using a mesh density function $h(R)=0.01\exp (0.5\sqrt{R})$ in DistMesh [33], does indeed improve the approximation of the energy density. However, this improvement comes at the cost of a significant increase in computational time from under 30 s for each of the simulations on the coarse mesh with 283 elements to around 5 min for the simulation on the locally refined mesh with 5859 elements. Hence, in order to obtain qualitatively similar results and to preserve the circular symmetry inherent in the problem, applying local mesh refinement leads to a more substantial increase in computational time than instead using a quadratic spatial basis approximation.

Fig. 5

Energy density (on a logarithmic scale) in a circular region of dissipative free-space with a point source at the origin (marked with a black star). Upper left simulation results with $N_s=0$ and $N_p=12$. Upper right simulation results with $N_s=2$ and $N_p=12$. Lower left simulation result with a locally refined mesh near the source point with $N_s=0$ and $N_p=12$. Lower right exact energy density plotted on the original (coarse) mesh. Parameter values: $\omega =100\pi$, $c=1$, $\rho ^f=1$ and $\mu =1$

Fig. 6

Energy density (on a logarithmic scale) inside a dissipative closed cavity for various orders of basis approximation. Parameter values: $\omega =100\pi$, $c=1$, $\rho ^f=1$, $\mu =0.5$ and source point (10, 5) (marked with a black star)

Fig. 7

Energy density (on a logarithmic scale) inside a dissipative closed cavity for two different order basis approximations. Parameter values: $\omega =100\pi$, $c=1$, $\rho ^f=1$, $\mu =1$ and source point (2.5, 7) (marked with a black star)

4.2 Closed Cavity Simulation

We perform numerical simulations on a closed polygonal cavity taken from Ref. [22]. In Fig. 6, we have placed a source point inside the cavity at (10, 5); whereas in Fig. 7, the source point is located on the boundary at (2.5, 7). The total number of mesh elements is 538 and we take $\omega =100\pi$, $c=1$ and $\rho ^f=1$ for all computations.

Figure 6 shows the predicted energy density for different approximation orders $N_s$ and $N_p$, and with dissipation rate $\mu =0.5$. The upper left plot of Fig. 6 shows the result with a piecewise constant approximation in space and $N_p=12$ in direction. Localized energy stripes near the source point are again evident, as they were in the free space propagation problem. The upper right plot of Fig. 6 shows the same calculation, but now with $N_s=2$. We notice that the circular symmetry around the source has been restored for only a modest increase in the computational time from around 30–45 s. The lower plots of Fig. 6 show the results of higher order simulations in both the space and the direction basis approximations with a maximum computational time of 6 min. The results are visually similar to those with $N_s=2$ and $N_p=12$. The sequence of plots in Fig. 6 demonstrate the convergence of DFM for increasing $N_s$ and the gain in accuracy that can be achieved by applying spatial approximations with $N_s\geqslant 2$ close to a source point.

The density distribution inside a cavity will depend on a number of factors including the geometry, the source point location and the damping. Figure 7 shows the result of moving the source point to the boundary and increasing the damping to $\mu =1$. Notice the clear difference between the energy density in the upper and lower parts of the cavity. The result in the left plot was computed using a piecewise constant approximation in space, and shows a sudden jump between the energy densities in the upper and lower regions. In the right plot we observe a smoother decay from the lower region to the upper region as a result of increasing the spatial approximation order. Therefore the right plot better captures the wave energy transmitted indirectly from the source to the upper part of the cavity due to reflections at the lower boundary. Hence, a higher order spatial approximation can restore symmetry and reduce spurious shadowing effects in both the near- and far-field of a point source excitation.

In all the numerical examples considered so far, the propagation has been through a homogeneous medium. We conclude this section by considering the case of a coupled two-cavity system with differing propagation speeds in each cavity.

Fig. 8

Simulations in a two-cavity system with propagation speeds $c_j=0.5$ if $x<0$ and $c_j=1$ otherwise, fluid densities $\rho ^f_{j}=k_j^{-2}$, and with hysteretic damping $\mu _j=\omega \eta /(2c_j)$, $\eta =0.01$ and $j=1,2,\ldots ,N$. Left energy density (on a logarithmic scale) in a coupled cavity system with source point $(-0.4,0.5)$ (marked with a black star), $\omega =100\pi$ and with $N_s=2$ and $N_p=24$. Right comparison of the energy ratio $\mathrm {P}_\text {left}/\mathrm {P}_\text {right}$ for different values of $N_s$ and $N_p$ in DFM compared with both DEA, and with 21 FEM simulations for the related wave problem at an equi-spaced range of frequencies within $\pm 5$ Hz of the ‘center’ frequencies used for the DEA and DFM computations

4.3 Coupled Two-Cavity System

Consider the two-cavity system shown in the left plot of Fig. 8 as considered in Ref. [7]. The source point is located at $(-0.4,0.5)$ and the energy density is computed for $N_s=2$ and $N_p=24$, with hysteretic damping applied via $\mu _j=\omega \eta /(2c_j)$, $\eta =0.01$, $\omega =100\pi$ and $\rho ^f_{j}=k_{j}^{-2}$ for $j=1,2,\ldots ,N$. We have taken $N=1184$ and set the propagation speeds to be $c_j=0.5$ when $\varOmega _j$ is within the left cavity ($x<0$) and $c_j=1$ otherwise. The total energy in each cavity is given by

$$\begin{aligned} \mathrm {P}_{\text {left}} = \int _{\varOmega |_{\left\{ x<0\right\} }} \rho _{\varOmega }(\mathbf {r}) \, \mathrm {d}\mathbf {r} \quad \text {and}\quad \mathrm {P}_{\text {right}} = \int _{\varOmega |_{\left\{ x>0\right\} }} \rho _{\varOmega }(\mathbf {r}) \, \mathrm {d}\mathbf {r}, \end{aligned}$$

where $\rho _{\varOmega }$ is the interior density (8) and the subscripts of $\mathrm {P}$ specify the cavity.

The right plot of Fig. 8 shows the ratio $\mathrm {P}_\text {left}/\mathrm {P}_\text {right}$ for five different frequency values. We compare the results of the DFM algorithm presented here for various orders of approximation with the DEA results given in [7]. The spectral convergence of the integrals arising in the DFM computation for this example is demonstrated in “Appendix”. The DEA results use a $6^{\mathrm {th}}$ order Chebyshev basis approximation along each edge of the polygonal cavity, and the same order of approximation globally in direction. The right plot of Fig. 8 also shows the results of 21 finite element method (FEM) simulations for the associated Helmholtz equation wave problem at an equi-spaced range of frequencies within $\pm 5$ Hz of the center frequencies used in the DEA and DFM computations. The FEM computations are performed using discontinuous Galerkin methods as reported in [7].

The results shown in the right plot of Fig. 8 demonstrate a good agreement between the DFM and the DEA results, particularly for lower damping (i.e. lower frequencies). We note that for finer triangle meshes than the one employed here (larger N), the results of the three DFM simulations become indistinguishable. For the mesh in the computations here, however, there is a considerable improvement when $N_s$ is increased from 0 to 2. The DFM results with $N_s=2$ and $N_p=12$ or $N_p=24$ agree well with the DEA result and lie towards the center of the range of FEM wave problem results. The computational times for these results were approximately 95 s per frequency when $N_s=0$, 150 s per frequency when $N_s=2$, $N_p=12$ and 550 s per frequency when $N_s=2$, $N_p=24$. Note that this compares favourably with the computational times for the DEA method reported in [7]; the 6${\text {th}}$ order DEA result shown in Fig. 8 took around 3000 s per frequency to compute.

4.4 Discussion: Computational Costs and Scaling

The computational costs of running a MATLAB based DFM code for the numerical experiments in this paper have been reported in the previous subsections. The times quoted are for non-optimised code running on a 3.4 GHz processor and without taking advantage of the fact that DFM is embarrassingly parallel with respect to the number of mesh elements. The cost of the algorithm as a whole is dominated by the cost of computing the entries of the matrix $\mathbf {B}$ and so scales with the number of non-zero entries of this matrix. As the problem size grows, the costs associated with solving the linear system become more significant, particularly in the case of low dissipation, but this still plays a relatively minor role for the examples presented here. An important advantage of the method over conventional numerical solvers for the Helmholtz equation when the wavenumber k becomes large is that the computational costs scale independently of k.

The approximate scaling of the algorithm as the number of degrees of freedom is increased is linear with respect to the number of mesh cells, since each triangular mesh element only transmits rays to a maximum of three other cells with which it shares a direct physical connection (as well as possible reflections into the same element). The scaling is less favourable with respect to the basis approximation order, and approximately scales as $O(N_s^2 N_p^2)$. However, we note that our results demonstrate a significant increase in accuracy for a moderate increase in the order of the spatial basis approximation from $N_s=0$ to $N_s=2$. In order to achieve the quadratic scaling of the algorithm with respect to the basis order in practice, then optimised integration routines such as the semi-analytic spectral methods described in Sect. 3.3 are crucial. Otherwise, the computational costs associated with the quadrature would also grow with the basis order in a sub-optimal manner. The results shown in Sect. 4 therefore demonstrate that DFM with p-refinement using semi-analytic integration and spectrally convergent quadrature is an efficient method for predicting energy distributions in complex structures.

5 Conclusions

A Discrete Flow Mapping algorithm with p-refinement in phase-space has been presented for approximating phase-space densities on triangulated domains. In particular, this is the first study incorporating orthogonal polynomial basis approximations of any specified order $N_s\ge 0$ in position space. The additional computational cost resulting from the higher order spatial approximations has been minimised by a careful evaluation of the propagation operator using semi-analytic integration methods for a Legendre polynomial basis, or full spectral quadrature in general. The numerical results presented verify the approach and demonstrate that for practical applications, a moderate increase in the order of the spatial approximation from $N_s=0$ to $N_s=2$ can be sufficient to improve the qualitative nature of the numerical solution.

Appendix: Spectrally Convergent Quadrature in Phase-Space

In this appendix we demonstrate the spectral convergence of the adaptive quadrature algorithm outlined in Sect. 3.3 using quadrature for both integrals. Consider two triangles connected along the line $x=0$ as shown in Fig. 9. These triangles are taken from the mesh of the coupled two-cavity system shown in Fig. 8. We give sample calculations for rays originating on the edge marked AB in Fig. 9, and hitting the edge BC. We fix the polynomial degrees in (21) to be $m=m'=2$ and $n=n'=12$. Then subdividing the integration domain for the direction integral as outlined in Sect. 3.1, we obtain the following four sub-intervals:

$$\begin{aligned} \mathrm {I}^{(1)} = (-0.6172,0.2095], \, \mathrm {I}^{(2)} = (0.2095,0.43],\, \mathrm {I}^{(3)} = (0.43,1.4772), \, \mathrm {I}^{(4)} = [1.4772,\pi /2). \end{aligned}$$

The first subdivision comes from the splitting (20), and the other subdivisions are necessary since the extremal values of the interval $\mathrm {I}^{(3)}$ correspond to branch points of the formula for $\theta _t$ obtained from Snell’s law (16). The integrals over $\mathrm {I}^{(1)}$, $\mathrm {I}^{(2)}$ and $\mathrm {I}^{(4)}$ are restricted to values of $\theta _j$ for which total internal reflection takes place ($\lambda _t=0$). For the integration over $\mathrm {I}^{(3)}$ we have both reflection and transmission, and the integration is performed with respect to $\theta _{t}^{+}$ as described in Sect. 3.1.

Fig. 9

Quadrature convergence study for mesh elements with different propagation speeds: $c_1=0.5$ and $c_2=1$. Vertex coordinates: $A=(-0.0456 ,0.4946)$, $B=(0,0.4622)$, $C=(0,0.5366)$ and $D=(0.0452,0.4988)$

Fig. 10

Spectral convergence of the quadrature rules in space and direction. Left integration over direction using four sub-divisions $\mathrm {I}^{(\kappa )}$, $\kappa =1,2,3,4$. Only the segment $\mathrm {I}^{(3)}$ includes non-trivial reflection/ transmission behaviour and a plot for each of the reflected and transmitted terms is given. Right integration over space at four fixed direction values. Parameter values: $m=m'=2$, $n,n'=12$, $c_1=0.5$, $c_2=1$, $\rho ^f_{1}=\rho ^f_{2}=1$ and $\mu _j=\omega \eta /(2c_j)$ with $\eta =0.01$ and $j=1,2$

The spectral convergence of the quadrature for the direction integration is shown in the left plot of Fig. 10. The plot shows the relative errors, where a numerical solution with $2^{12}+1$ quadrature nodes is used in the place of an exact solution. Notice that the integrals with total internal reflection (over $\mathrm {I}^{(1)}$, $\mathrm {I}^{(2)}$, and $\mathrm {I}^{(4)}$) converge faster than the integrals with a non-constant transmission probability function $\lambda _t$. In the right plot of Fig. 10, we also demonstrate the spectral convergence of the spatial quadrature. We consider the spatial integrals with four fixed $\theta _j$ values. The choices of $\theta _j$ shown in the plot correspond to the central value for each sub-interval. Note that the relative error for these integrals has already reached single precision with just 8 quadrature nodes, and double precision with 16 quadrature nodes. This rapid convergence means that the computational cost depends most crucially on the efficiency of the direction coordinate integration, since this specifies how many spatial integrals must be evaluated. The spatial integration efficiency could only practically be improved by the use of exact spatial integration formulae, such as (23), or (25) for $N_s>0$.