Tent pitching and Trefftz-DG method for the acoustic wave equation

TL;DR

This paper introduces a space-time Trefftz discontinuous Galerkin method for efficiently solving the acoustic wave equation on tent pitched meshes, enabling local and parallel computations in heterogeneous media.

Contribution

It develops a semi-explicit DG Trefftz method on tent pitched meshes, with implementation details and parallelization strategies for heterogeneous media.

Findings

Efficient local solution of acoustic wave equation using tent pitched meshes.

Parallel implementation in NGSolve enhances computational speed.

Method effectively handles heterogeneous media.

Abstract

We present a space-time Trefftz discontinuous Galerkin method for approximating the acoustic wave equation semi-explicitly on tent pitched meshes. DG Trefftz methods use discontinuous test and trial functions, which solve the wave equation locally. Tent pitched meshes allow to solve the equation elementwise, allowing locally optimal advances in time. The method is implemented in NGSolve, solving the space-time elements in parallel, whenever possible. Insights into the implementational details are given, including the case of propagation in heterogeneous media.

Tables

Table 1. Table 1 : Convergence rates and run time comparison for a singular solution on the L-shapes domain, comparing uniform meshing and meshes refined towards the singularity.

mesh	$h_{\max }$	total #dofs	L2-error	dof-rate	runtime [s]
uniform	0.07	$\mathrm{320\hspace{0.17em}112}$	$0.0178$	-	0.6234
	0.05	$\mathrm{875\hspace{0.17em}696}$	$0.0117$	1.2638	1.5916
	0.03	$\mathrm{3\hspace{0.17em}894\hspace{0.17em}752}$	$0.0067$	1.1022	7.9157
	0.01	$\mathrm{102\hspace{0.17em}320\hspace{0.17em}048}$	$0.0021$	1.0662	255.3233
non-uniform	0.12	$\mathrm{1\hspace{0.17em}138\hspace{0.17em}256}$	$0.0218$	-	3.276
	0.10	$\mathrm{2\hspace{0.17em}003\hspace{0.17em}536}$	$0.0083$	5.1308	4.6809
	0.08	$\mathrm{3\hspace{0.17em}832\hspace{0.17em}624}$	$0.0030$	4.7041	8.0069
	0.06	$\mathrm{9\hspace{0.17em}777\hspace{0.17em}536}$	$0.0008$	4.2104	23.4588

Equations

$$\{\begin{array}{ll}\nabla \cdot \mathbf{\sigma }+c^{-2}\frac{\partial v}{\partial t}=0& \text{ in }Q,\\ \nabla v+\frac{\partial \mathbf{\sigma }}{\partial t}=0& \text{ in }Q,\\ v(\cdot ,0)=v_0,\mathbf{\sigma }(\cdot ,0)={\mathbf{\sigma }}_0& \text{ on }\Omega ,\\ v=g_D& \text{ on }{\Gamma }_D\times [0,T],\\ {\mathbf{n}}_{\Omega }^x\cdot \mathbf{\sigma }=g_N& \text{ on }{\Gamma }_N\times [0,T],\end{array}$$

$$\{\begin{array}{ll}-\Delta U+c^{-2}\frac{{\partial }^{2}U}{\partial t^2}=0& \text{ in }Q,\\ \frac{\partial U}{\partial t}(\cdot ,0)=v_0,U(\cdot ,0)=U_0& \text{ on }\Omega ,\\ \frac{\partial U}{\partial t}=g_D& \text{ on }{\Gamma }_D\times [0,T],\\ -{\mathbf{n}}_{\Omega }^x\cdot \nabla U=g_N& \text{ on }{\Gamma }_N\times [0,T].\end{array}$$

$$c|{\mathbf{n}}_F^x|<n_F^t\text{ and we call the face space-like, or}$$

$$n_F^t=0\text{ and we call the face time-like.}$$

$${\mathcal{F}}_h$$

$${\mathcal{F}}_h^{\text{space}}$$

$${\mathcal{F}}_h^{\text{time}}$$

$${\mathcal{F}}_h^0$$

$${\mathcal{F}}_h^D$$

$$\displaystyle\bm{T}(K):=\big{\{}(w,\bm{\tau})\in L^{2}(K)^{1+n}\text{ s.t. }{\left.\kern-1.2pt\bm{\tau}\right|_{\partial K}}\in L^{2}(\partial K)^{n},\ \frac{\partial{w}}{\partial t},\nabla\cdot\bm{\tau}\in L^{2}(K),$$

$$\displaystyle\frac{\partial{\bm{\tau}}}{\partial t},\nabla w\in L^{2}(K)^{1+n},\ \nabla\cdot\bm{\tau}+c^{-2}\frac{\partial w}{\partial t}=0,\ \nabla w+\frac{\partial\bm{\tau}}{\partial t}=\bm{0}\ \big{\}}$$

$$\displaystyle\bm{T}(\mathcal{T}_{h}):=\Big{\{}(w,\bm{\tau})\in L^{2}(Q)^{1+n},\ \text{ s.t. }({\left.\kern-1.2ptw\right|_{K}},{\left.\kern-1.2pt\bm{\tau}\right|_{K}})\in\bm{T}(K)\ \forall K\in\mathcal{T}_{h}\Big{\}},$$

$$\begin{array}{rl}& -{\int }_{K}v\left(\nabla \cdot \mathbf{\tau }+c^{-2}\frac{\partial w}{\partial t}\right)+\mathbf{\sigma }\cdot \left(\frac{\partial \mathbf{\tau }}{\partial t}+\nabla w\right)dV\\ & +{\int }_{\partial K}v\left(\mathbf{\tau }\cdot {\mathbf{n}}_K^x+c^{-2}w{n}_K^t\right)+\mathbf{\sigma }\cdot \left(w\cdot {\mathbf{n}}_K^x+\mathbf{\tau }{n}_K^t\right)dS=0.\end{array}$$

$$\begin{array}{rl}& {\int }_{\partial K}{\hat{v}}_{hp}\left(\mathbf{\tau }\cdot {\mathbf{n}}_K^x+c^{-2}w{n}_K^t\right)+{\hat{\mathbf{\sigma }}}_{hp}\cdot \left(w\cdot {\mathbf{n}}_K^x+\mathbf{\tau }{n}_K^t\right)dS=0.\end{array}$$

$$\{\{w\}\}$$

$$[[w]{]}_t$$

$$[[w]{]}_N$$

$${\hat{v}}_{hp}=\{\begin{array}{l}{v}_{hp}^{-}\\ v_0\\ \{\{v_{hp}\}\}+\beta [[{\mathbf{\sigma }}_{hp}]{]}_N\\ v_{hp}\\ g_D\\ v_{hp}+\beta ({\mathbf{\sigma }}_{hp}\cdot {\mathbf{n}}_{\Omega }^x-g_N)\end{array}{\hat{\mathbf{\sigma }}}_{hp}=\{\begin{array}{ll}{\mathbf{\sigma }}_{hp}^{-}& \text{on }{\mathcal{F}}_h^{\text{space}},\\ {\mathbf{\sigma }}_0& \text{on }{\mathcal{F}}_h^0,\\ \{\{{\mathbf{\sigma }}_{hp}\}\}+\alpha [[v_{hp}]{]}_N& \text{on }{\mathcal{F}}_h^{\text{time}},\\ {\mathbf{\sigma }}_{hp}& \text{on }{\mathcal{F}}_h^T,\\ {\mathbf{\sigma }}_{hp}-\alpha (v_{hp}-g_D){\mathbf{n}}_{\Omega }^x& \text{on }{\mathcal{F}}_h^D,\\ g_N{\mathbf{n}}_{\Omega }^x& \text{on }{\mathcal{F}}_h^N,\end{array}$$

$$\begin{array}{rl}& \text{find }(v_{hp},{\mathbf{\sigma }}_{hp})\in {\mathbf{V}}_p({\mathcal{T}}_h)\text{s.t.}\\ & \mathcal{A}(v_{hp},{\mathbf{\sigma }}_{hp};w,\mathbf{\tau })=\ell (w,\mathbf{\tau })\forall (w,\mathbf{\tau })\in {\mathbf{V}}_p({\mathcal{T}}_h),\end{array}$$

$$\mathcal{A}(v_{hp},{\mathbf{\sigma }}_{hp};w,\mathbf{\tau })$$

$$+{\int }_{{\mathcal{F}}_h^{\text{time}}}\left(\{\{v_{hp}\}\}[[\mathbf{\tau }]{]}_N+\{\{{\mathbf{\sigma }}_{hp}\}\}\cdot [[w]{]}_N+\alpha [[v_{hp}]{]}_N\cdot [[w]{]}_N+\beta [[{\mathbf{\sigma }}_{hp}]{]}_N[[\mathbf{\tau }]{]}_N\right)dS$$

$$+{\int }_{{\mathcal{F}}_h^T}{c}^{-2}{v}_{hp}w+{\mathbf{\sigma }}_{hp}\cdot \mathbf{\tau }dS+{\int }_{{\mathcal{F}}_h^D}(\mathbf{\sigma }\cdot {\mathbf{n}}_{\Omega }^x+\alpha v_{hp})wdS$$

$$+{\int }_{{\mathcal{F}}_h^N}{v}_{hp}(\mathbf{\tau }\cdot {\mathbf{n}}_{\Omega }^x)+\beta (\mathbf{\sigma }\cdot {\mathbf{n}}_{\Omega }^x)(\mathbf{\tau }\cdot {\mathbf{n}}_{\Omega }^x)dS$$

$$\ell (w,\mathbf{\tau })$$

$${\mathbb{T}}^p(K)$$

$$\dim {\mathbb{T}}^p(K)=(n+1)\left(\begin{array}{c}p+n\\ n\end{array}\right)={\mathcal{O}}_{p\to \infty }(p^n),$$

$${\mathbb{U}}^p(K):=\left\{U\in {\mathbb{P}}^p(K):-\Delta U+\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}U=0\right\}.$$

$$U(\mathbf{x},t)=\sum _{\begin{array}{c}\mathbf{\alpha }\in {\mathbb{N}}_0^n,k\in {\mathbb{N}}_0,\\ |\mathbf{\alpha }|+k\le p\end{array}}{a}_{k,\mathbf{\alpha }}{\mathbf{x}}^{\mathbf{\alpha }}{t}^k.$$

$$a_{k,\mathbf{\alpha }}=\frac{c^2}{k(k-1)}\sum _{m=1}^n({\alpha }_m+1)({\alpha }_m+2)a_{k-2,\mathbf{\alpha }+2{\mathbf{e}}_m}$$

$$\dim {\mathbb{U}}^p(K)=\left(\begin{array}{c}p+n\\ n\end{array}\right)+\left(\begin{array}{c}p-1+n\\ n\end{array}\right)=\frac{2p+n}{p}\left(\begin{array}{c}p-1+n\\ n\end{array}\right).$$

$${\mathbb{U}}^{p+1}(K)=\text{span}\{b_j,j\in \mathcal{I}\}\text{ by setting }{\mathbb{W}}^p(K)=\text{span}\{(\frac{\partial b_j}{\partial t},-\nabla b_j),j\in \mathcal{I}\}.$$

$$\dim {\mathbb{W}}^p(K)=\dim {\mathbb{U}}^{p+1}(K)-1=\frac{2p+n+2}{p+1}\left(\begin{array}{c}n+p\\ n\end{array}\right)$$

$$h_K:=\underset{(\mathbf{x},t),(\mathbf{y},s)\in K}{\sup }(|\mathbf{x}-\mathbf{y}{|}^2+c^2|t-s{|}^2{)}^{1/2}$$

$$(\mathbf{x},t)=(h_K\hat{\mathbf{x}},h_K{c}^{-1}\hat{t})$$

$$v(\mathbf{x},t)=c\hat{v}(\hat{\mathbf{x}},\hat{t}),\mathbf{\sigma }(\mathbf{x},t)=\hat{\mathbf{\sigma }}(\hat{\mathbf{x}},\hat{t}).$$

Full text

††journal: Computers & Mathematics with Applications

\PreventPackageFromLoading

natbib

\useunder\ul

1 Introduction

Standard finite element methods approximate the solution of a given partial differential equation (PDE) by piecewise polynomial functions. A classic approach to the discretization of time dependent PDEs is to use finite element methods to discretize space and then use time stepping schemes to advance in time. We consider a different approach based on finite element approximation simultaneously in space and time. In this paper, we present original numerical results for the space-time Trefftz discontinuous Galerkin (DG) method studied in [20], confirming the known, and conjectured, properties of the method numerically. This is the first implementation combining a Trefftz-DG method with tent pitched meshes, highlighting the importance of both techniques in order to produce an efficient method. We point out that we can solve the problem in $n+1$ dimensions, $n\in\{1,2,3\}$, without the need for $n+1$ dimensional elements, since the Trefftz method only requires the computation of integrals at interelement boundaries. Furthermore, we introduce a way of recovering the solution of the second order system, and address in detail the issue of inhomogeneous materials.

The space-time approach requires us to mesh the full space-time domain. Furthermore, the use of approximation spaces based on piecewise ”total degree” polynomials in both space and time, leads to a higher number of degrees of freedom. On the upside, $hp$-refinement is made possible in space-time, allowing for straightforward higher order approximation. Additionally, the space-time domain mesh is not forced to be a product mesh, as it is for time stepping schemes. Instead, we are allowed to use unstructured meshes. This gives us the possibility to devise suitable mesh design strategies in order to circumvent the CFL-condition, which usually limits the global time-step size by the size of the smallest spatial element in explicit time stepping schemes. This will allow to advance in time elementwise and in parallel.

Space-time finite element methods for linear wave propagation go as far back as [15], and have been used with DG methods e.g. in [21, 9, 18, 5]. DG methods are based on discontinuous piecewise polynomial functions, and so-called numerical fluxes which impose continuity constraints at mesh inter-element boundaries.

On the one hand, ideas on combining DG methods with so called tent pitched meshes can be found in e.g. [22, 21, 9, 19, 29]. A way of constructing these meshes can, for instance, be found in [27, 8, 12, 28]. Tent pitching techniques generate a space-time mesh, which complies with the causality properties of the hyperbolic PDE. The resulting mesh consists of tent shaped objects, each advancing locally optimally in time, with the PDE being explicitly solvable in each of them. Though the tent pitching strategy pairs well with DG methods, also other methods are applicable in combination with tent pitched meshes. In [12, 28, 10], schemes for the semi-discretization of different hyperbolic equations on tents are presented, which map tents to a domain where space and time are separable. Similar to the Trefftz-DG method, these schemes are able to solve in 3+1 dimensions, without building four dimensional elements. Friedrichs theory is used in [11] in order to derive a conforming methods, and to prove its convergence properties. We point out that tent pitching is not the only way to deal with the time step restriction of locally refined meshes. A stabilization for a conforming space-time finite element method on Cartesian (in time) meshes is presented in [25]. Classical time-stepping schemes can still be applied successfully by splitting the domain into a coarse-mesh and a fine-mesh region, then explicit time stepping in the coarse-mesh region is combined with local implicit or explicit time stepping in the fine-mesh region. A fully explicit scheme can be found in [14, 13].

On the other hand, Trefftz methods, originating from [26], incorporate properties of the PDE into the test and trial spaces. This is done by choosing them as (local) solutions of the targeted differential equation. The use of Trefftz spaces allows us to reduce the number of degrees of freedom, as compared to the total degree polynomial spaces, however keeping the same accuracy. Work on Trefftz-DG methods for different wave propagation problems includes [16, 3, 4, 7, 17, 20, 6]. In [3], a Trefftz-DG method in space-time for the second order wave equation is presented, proving $h$-convergence in 1, 2 and 3 space dimensions, as well as $hp$-convergence, along with exponential convergence for analytic solutions in 1 space dimension. In [4], Trefftz-DG is applied to the coupled elasto-acoustic system, and well-posedness of the problem, as well as error estimates in mesh-dependent norms, are shown. Both [3, 4] are formulated for meshes with tensor product structure in time. In [20] a Trefftz-DG method for the acoustic wave equation in first order formulation and in arbitrary space dimension is presented. The formulation works for tensor product (in time) meshes, as well as for tent pitched meshes. Well-posedness and optimal $h$-convergence are proven, but no numerical results are shown.

In this paper, we focus on the combination of Trefftz-DG formulation and tent pitched meshes, and on its efficient implementation. As Trefftz-DG formulations only contain interelement terms, they pose a natural choice to evolve the solution from the bottom to the top of tent elements. As mentioned, previous implementations on Trefftz-DG methods were based on tensor product structure in time [3, 4, 17].

This work proceeds as follows. First, we introduce the Trefftz-DG method in Section 2, starting by stating the model problem, defining the Trefftz spaces and finishing the section by formulating the method, as it was introduced in [20]. We continue in Section 3 by reviewing different strategies of discretizing the Trefftz spaces. In Section 4, we discuss some numerical details on how to evolve the solution elementwise on a tent pitched mesh, and in Section 5, we show a way to recover the second order solution from the first order formulation. Finally, we present numerical results, which were obtained by implementation of the method in NGSolve [23, 24], in Section 6.

2 The Trefftz-DG method for the acoustic wave equation

2.1 The acoustic wave equation

Let $Q=\Omega\times(0,T)$ be a space-time domain in $\mathbb{R}^{n+1},\ n\in\{1,2,3\}$, where $\Omega\subset\mathbb{R}^{n}$ is a Lipschitz bounded domain with outward unit normal vector $\bm{n}^{\bm{x}}_{\Omega}\in\mathbb{R}^{n}$. The boundary of $\Omega$ is divided into two parts $\Gamma_{D}$ and $\Gamma_{N}$, corresponding to Dirichlet and Neumann boundary conditions, respectively, such that they have disjoint interiors, the union of their closures gives the boundary of $\Omega$, and one of them can be empty. We consider the acoustic wave equation in first order formulation, given by

[TABLE]

where we assume that the wavespeed $c>0$ is piecewise constant on $\Omega$.

If the initial condition $\bm{\sigma}_{0}$ it the gradient of a scalar field $U_{0}$, i.e. $\bm{\sigma}_{0}=-\nabla U_{0}$, then the first order system is equivalent to the second order system obtained by setting $v=\frac{\partial{U}}{\partial t}$ and $\bm{\sigma}=-\nabla U$:

[TABLE]

The Laplacian $\Delta$, gradient $\nabla$ and divergence $\nabla\cdot$ are considered with respect to the space variable $\bm{x}$ only.

2.2 Space-time meshes

The mesh $\mathcal{T}_{h}(Q)$ of the space-time domain $Q$ is assumed to consist of non-overlapping Lipschitz polytopes, where $h=\max_{K\in\mathcal{T}_{h}(Q)}h_{K}$, with $h_{K}$ being the anisotropic diameter defined in (3.2). For each mesh face $F=\partial K_{1}\cap\partial K_{2}$, for $K_{1},K_{2}\in\mathcal{T}_{h}(Q)$, we assume that it either lies below the characteristic speed $1/c$, or is vertical (parallel to the time axis). In more rigorous terms: Let $(\bm{n}_{F}^{x},n_{F}^{t})$ be the normal vector to $F$ with $n_{F}^{t}\geq 0$, then either

[TABLE]

Notice, however, that no CFL-condition or any other time step size restriction is imposed on the time-like faces.

A mesh with space-like faces only is called a tent pitched mesh. In the numerical experiments presented below, we generate tent pitched meshes by using the algorithm presented in [12] (see also [28]). The mesh is buit by progressively advancing in time, stacking tent-shaped objects on top or each other, each of them union of $(n+1)$-simplices. The main idea is that the tent height is chosen such that the differential equation is explicitly solvable in each tent. Therefore, the local maximal time advance at a spatial point has to respect the causality constraint, which corresponds to a local CFL-condition. This allows us to advance the solution tent by tent, not necessarily having to solve a global system. For independent tents, i.e. tents that are not on top of each other, the computations can be done in parallel. We remark that, in the numerical tests, we observe no stability issues with tents pitched very close to the limit of the causality condition.

The set of all faces of our mesh,

[TABLE]

can be separated into two sets of internal faces

[TABLE]

and sets of faces on the boundary of $Q$, split up according to their types of initial and boundary conditions

[TABLE]

This classification of the faces is represented in Figure 1.

2.3 Trefftz spaces

By definition, Trefftz functions in the kernel of the considered differential operator. For the first order wave equation, we define the local and global Trefftz space as

[TABLE]

[TABLE]

respectively. Note that, by assuming that the solution is in $\bm{T}(\mathcal{T}_{h})$, we require additional smoothness on the solution, as in general we only have that $\frac{\partial{v}}{\partial t}+\nabla\cdot\bm{\sigma}\in L^{2}(K)$, for all $K\in\mathcal{T}_{h}$.

We derive the Trefftz-DG method for any choice of discrete test and trial space with a Trefftz property, which we denote by $\bm{V}_{p}(\mathcal{T}_{h})$. A possible choice for a polynomial $\bm{V}_{p}(\mathcal{T}_{h})$ is given in Section 3 below.

2.4 The Trefftz-DG method

Following [20], we derive the Trefftz-DG method for the IBVP in (2.1). The method is derived from a local weak formulation, obtained by multiplying the two equations in (2.1) by test and trial functions $w$ and $\bm{\tau}$, respectively, and integrating by parts on each element $K$ of the mesh $\mathcal{T}_{h}(Q)$. Then, adding the two equations gives

[TABLE]

By choosing Trefftz test functions $(w,\bm{\tau})\in\bm{V}_{p}(K)$, the volume integrals over $K$ vanishes. We are left with:

[TABLE]

Typical for DG methods, the continuity of the numeric solution on inter-element boundaries is enforced within the bilinear form of the method. To this end, the trace of the solution $(v,\bm{\sigma})$ in the boundary integral has been replaced by the numeric fluxes $(\hat{v}_{hp},\hat{\bm{\sigma}}_{hp})$, which we define below.

To do so, we need to introduce some standard DG notation. For a face $F=\partial K^{+}\cap\partial K^{-}$ shared by two elements $K^{+},K^{-}\in\mathcal{T}_{h}(Q)$, we define the average $\{\!\!\{\cdot\}\!\!\}$, jumps in space $[\![\cdot]\!]_{N}$, and jumps in time $[\![\cdot]\!]_{t}$, for scalar- and vector-valued functions as follows

[TABLE]

where $n_{K}=(\bm{n}_{K}^{x},n_{K}^{t})$ is the unit outer normal vector at $\partial K$ split into its space and time components.

Across time-like faces, the information is passed by using centered fluxes with jump penalization, whereas, across space-like faces, the information is passed upward in time, resembling an up-wind scheme. More precisely, the fluxes on the inter-element faces are chosen as

[TABLE]

where $\alpha$ and $\beta$ are penalty parameters, which will be chosen constant (notice that they are needed on time-like and Dirichlet faces only). By $w^{+}$ and $w^{-}$ we denote the trace of the function $w$ on space-like faces from the adjacent element at higher and lower times, respectively.

Finally, we plug the definition of the fluxes into (2.4) and sum over all elements $K\in\mathcal{T}_{h}(Q)$. Then the Trefftz-DG method for the wave equation reads:

[TABLE]

with

[TABLE]

and

[TABLE]

On a tent pitched mesh, as the one in Figure 1, the method is semi-explicit, meaning that the solution on each tent only depends on the tents below, allowing to solve each tent explicitly, and tents independent from each other in parallel; details are given in Section 4.1 below. The situation where also vertical faces are present, is needed, for instance, in the case of piecewise constant wavespeed, is discussed in Section 4.2 below. Note that the method only includes integrals over element boundaries, thus only quadrature on $n$ dimensional simplices is needed.

3 Choice of discrete Trefftz spaces

So far, we have not specified what discretization of the Trefftz space $\bm{V}_{p}(K)\subset\bm{T}(K)$ to use. We introduce the straightforward choice, given by all polynomials in space-time that fulfill the first order wave equation. For an element $K\subset\mathbb{R}^{n+1}$ in the mesh $\mathcal{T}_{h}(Q)$, we define the local polynomial Trefftz space as

[TABLE]

where we denote by $\mathbb{P}^{p}(K)$ the space of polynomials on $K$ of degree $\leq p$. In general, it is possible to choose different polynomial degrees $p$ in different elements. Here, we choose a uniform $p$, as this is consistent with the numerical examples below. The global Trefftz-DG space on the whole mesh is then given by $\mathbb{T}^{p}(\mathcal{T}_{h}):=\prod_{K\in\mathcal{T}_{h}}\mathbb{T}^{p}(K).$ The dimension of the elemental Trefftz space is given by

[TABLE]

where we recall that $\left(\begin{smallmatrix}a\\ b\end{smallmatrix}\right)=\frac{a!}{b!(a-b)!}$ for $b\leq a\in\mathbb{N}_{0}$. Notice that, for the total degree polynomial space, one has $\dim(\mathbb{P}^{p}(\mathbb{R}^{n+1})^{n+1})=\mathcal{O}_{p\rightarrow\infty}(p^{n+1})$.

Let us now assume that the first order problem is derived from a second order problem. Then it is natural to derive the vector valued Trefftz space for the first order problem from a scalar Trefftz space for the second order problem. We now detail this approach as it is the one we use for the numerical results presented in Section 6. Let us start by defining the polynomial Trefftz space for the second order problem:

[TABLE]

We are able to construct a basis for this space using the recursion formula introduced in [20, Remark 13]. We recall it here, for completeness. We need some multi-index notation: for $\bm{\alpha}\in\mathbb{N}^{n}_{0}$ we denote $|\bm{\alpha}|=\alpha_{1}+\dots+\alpha_{n}$ and $\bm{x}^{\bm{\alpha}}=x_{1}^{\alpha_{1}}\dots x_{n}^{\alpha_{n}}$. Furthermore, let $\bm{e}_{m}:=(0,\dots,0,1,0,\dots,0)\in\mathbb{N}^{n}_{0}$ with 1 in the $m$-th entry. Consider a space-time polynomial

[TABLE]

We want to compute the coefficients $a_{k,\bm{\alpha}}$ such that the polynomial is Trefftz. This is done by inserting the polynomial into the second order wave equation and collecting terms of equal power to find that

[TABLE]

has to hold for the polynomial to be Trefftz. To start the recursion, we need to choose polynomial bases (in the space variables only) for $k=0$ and $k=1$, respectively. More precisely, we start by choosing polynomial basis functions $\{\tilde{b}_{1},\dots,\tilde{b}_{\left(\begin{smallmatrix}p+n\\ n\end{smallmatrix}\right)}\}$ for the space $\mathbb{P}^{p}(\mathbb{R}^{n})$ and $\{\hat{b}_{1},\dots,\hat{b}_{\left(\begin{smallmatrix}p-1+n\\ n\end{smallmatrix}\right)}\}$ for $\mathbb{P}^{p-1}(\mathbb{R}^{n})$. Then we can introduce a basis for $\mathbb{U}^{p}(K)$ such that either $U(\cdot,0)=\tilde{b}_{\ell}$ and $\frac{\partial{U}}{\partial t}(\cdot,0)=0$, or $U(\cdot,0)=0$ and $\frac{\partial{U}}{\partial t}(\cdot,0)=\hat{b}_{\ell}$ for some $\ell$. Hence, we can construct the basis for $\mathbb{U}^{p}(K)$ out of two sets of polynomial basis functions of $\mathbb{P}^{p}(K)$ and $\mathbb{P}^{p-1}(K)$. This lets us determine the dimension as

[TABLE]

Then, a Trefftz space for the first order system can be derived from

[TABLE]

We have that

[TABLE]

and furthermore $\mathbb{W}^{p}(K)\subset\mathbb{T}^{p}(K)$. A recursion formula, similar to (3.1), can also be derived for $\mathbb{T}^{p}(K)$, however the numerical results in Section 6 are centered around $\mathbb{W}^{p}(K)$.

Remark 3.1**.**

It is sufficient to compute the coefficients only once for $c=1$ and then fix the wavenumber by a coordinate transform. Furthermore, for numerical stability, it is convenient to shift the basis functions to the center of the element and scale them by its anisotropic diameter, which is defined by

[TABLE]

for a mesh element $K$. For reference coordinates $(\hat{\bm{x}},\hat{t})$, the coordinate transform given by

[TABLE]

transforms the Trefftz basis $\hat{U}(\hat{\bm{x}},\hat{t})$ of wavespeed 1 to Trefftz basis functions $\hat{U}(x,t)$ of arbitrary wavespeed $c$. In the case of Trefftz functions for the first order system $(\hat{v},\hat{\bm{\sigma}})$, we need to choose

[TABLE]

4 Evolution within a tent

The tent pitched mesh allows us to solve local tents explicitly. This is due to the fact that the slope of the mesh faces is below the characteristic speed $1/c$, thus the local solution on a tent can be computed once the solution on its inflow boundary is known. In Section 4.1, we discuss how to evolve the solution within a tent with constant wavespeed inside the tent itself. The case where the wavespeed changes within a tent is considered in Section 4.2. Notice that, in the constant wavespeed case, tents coincide with mesh elements, while in the latter case tents on the interface contain more than one mesh element.

4.1 Constant wavespeed

Let us denote the bottom and top faces of the tent by ${T^{\text{bot}}_{h}\subset(\mathcal{F}_{h}^{\text{space}}\cup\mathcal{F}_{h}^{0})}$ and ${T^{\text{top}}_{h}\subset(\mathcal{F}_{h}^{\text{space}}\cup\mathcal{F}_{h}^{T})}$, respectively. Furthermore, tent faces on the boundary are denoted by $T^{D}_{h}\subset\mathcal{F}^{D}_{h}$ for Dirichlet and $T^{N}_{h}\subset\mathcal{F}^{N}_{h}$ for Neumann boundaries.

Since the solution is explicit on each tent, we only need to solve a local system of size ${\dim(\bm{V}^{p}(K))\times\dim(\bm{V}^{p}(K))}$. The system is derived from (LABEL:eq:TDG) and is given by the following equation

[TABLE]

where, in the case $T_{h}^{\text{bot}}\subset\mathcal{F}_{h}^{0}$, $(v_{\text{bot}},\bm{\sigma}_{\text{bot}})=(v_{0},\bm{\sigma}_{0})$, and in the case $T_{h}^{\text{bot}}\subset\mathcal{F}_{h}^{space}$, $(v_{\text{bot}},\bm{\sigma}_{\text{bot}})$ on a given face is the previously computed solution in the tent sharing that face in lower time.

For the numerical integration, we only need an integration rule for $n$-simplices, in order to integrate over the boundary of the tent. We can define an integration rule on the spatial mesh once, which we can then map to the faces of the tent. This idea is visualized in Figure 2. After solving on the tent, we need to evaluate $(v_{hp},\bm{\sigma}_{hp})$ in the integration points on $T^{\text{top}}_{h}$, and store these values for the next tent. On each spatial integration point, we only need to store the most recent results, leading to a total storage of: ${(\text{total number of integration points})\cdot(n+1)}$.

4.2 Piecewise constant wavespeed

Recall that we assume that the wavespeed is constant in time and piecewise constant in space. In this case, we always consider initial spatial meshes that are aligned with the discontinuities of the wavespeed. To treat such a jump within a space-time tent, we need to incorporate the jump terms from our DG formulation (LABEL:eq:TDG). This involves integrating on the time-like inter-element boundary contained inside the tent, denoted by $T^{\text{time}}_{h}\subset\mathcal{F}^{\text{time}}_{h}$. According to (LABEL:eq:TDG), one has to add to the left-hand side of (4.1) the term

[TABLE]

Since the tent now includes two mesh elements, the system matrix is now of size ${2\dim(\bm{V}^{p}(K))\times 2\dim(\bm{V}^{p}(K))}$. The extension to interfaces between more than two materials follows.

5 Recovery of the solution of the second order equation

In the case where the problem comes from a second order formulation we can substitute $v=\frac{\partial{U}}{\partial t}$ and $\bm{\sigma}=-\nabla U$ to write the method in terms of test and trial functions from $\mathbb{U}_{p+1}(K)$. Then the method (LABEL:eq:TDG) reads:

[TABLE]

with

[TABLE]

The constant basis function does not contribute to formulation (LABEL:eq:TDG2), as only derivatives of the unknown $U_{hp}$ are present. Thus, this formulation produces the same results as (LABEL:eq:TDG) with $\bm{V}_{p}(\mathcal{T}_{h})=\mathbb{W}_{p}(\mathcal{T}_{h})$. In order to fix the constants and recover the solution to the second order wave equation, we modify the original formulation by adding the additional terms

[TABLE]

to the bilinear form $\hat{\mathcal{A}}(U_{hp};V)$, and

[TABLE]

to the right hand side $\hat{\ell}(V)$, where and $U_{0}(x)=U(x,0)$. Note that these terms preserve the consistency of the formulation.

Therefore, when evolving the solution inside a single tent, we need to add $\int_{T^{\text{top}}_{h}}U_{hp}V$ and $\int_{T^{\text{bot}}_{h}}U_{\text{bot}}V$ to the left- and right-hand side, respectively, of the formulation discussed in Section 4.

6 Numerical results

In this section we present numerical test results in one, two, and three spatial dimensions. The Trefftz-DG method was implemented in NGSolve [23, 24]. If not otherwise stated, we use the following settings for the numerical examples. We consider the problem (2.1) with initial and Dirichlet boundary conditions such that the analytical solution is $(v,\bm{\sigma})=(\frac{\partial{U}}{\partial t},-\nabla U)$, where $U$ is the standing wave

[TABLE]

given here in 3+1 dimensions, and set the wavespeed $c=1$. An example is plotted in 1+1 dimensions in Figure 3. The penalty parameters are chosen as $\alpha=\beta=0.5$. We measure the error

[TABLE]

at final time $T$, which we choose at $T=1$. The tent pitched meshes are produced by the algorithm presented in [12]. In this algorithm, the height of the tents is limited by the slope of the edges, and not by the slope of the faces. All timings were performed on a server with two Intel(R) Xeon(R) CPU E5-2687W v4, with 12 cores each.

6.1 Approximation properties of Trefftz spaces

Independently of their combination with tent pitching, Trefftz polynomial spaces possess good approximation properties for wave solutions. In [20, Section 6] , $h$-version approximation estimates for wave solutions in Trefftz polynomial spaces were proven. The derivation of $p$-version approximation estimates in higher dimensions is, to the best of our knowledge, still open. Here, we compare the Trefftz space to a full discontinuous polynomial space. Recall that, thanks to the Trefftz property of the test functions, we were able to cancel the volume integral in the weak formulation (2.4). This does not hold for the full polynomial space. Therefore, we need to add the volume term to the left-hand side of the formulation (LABEL:eq:TDG), giving the new left-hand side:

[TABLE]

We now need to solve $\mathcal{\tilde{A}}(v_{hp},\bm{\sigma}_{hp};w,\bm{\tau})=\ell(w,\bm{\tau}),\ \forall(w,\bm{\tau})\in\mathbb{P}^{p}(\mathcal{T}_{h})^{n+1}$ for $(v_{hp},\bm{\sigma}_{hp})\in\mathbb{P}^{p}(\mathcal{T}_{h})^{n+1}$. Notice that, as opposed to the Trefftz-DG method (LABEL:eq:TDG) the (full) DG method requires the computation of integrals also in space-time volumes.

For the numerical results we have taken the unit square domain, uniformly meshed into space-time squares, and used the Trefftz space $\mathbb{W}^{p}$ or the space $\mathbb{Q}^{p}$ of polynomials with maximal order $p$ in each component. In the latter case the number of degrees of freedom per element is hence given by $2(p+1)^{2}$.

The results are shown in Figure 4. Both choices exhibit similar, exponential, convergence speed in terms of polynomial degree, although the Trefftz space is only a subset of the polynomials of maximal degree equal to $p$. The benefits of the Trefftz space becomes clear when comparing errors versus number of degrees of freedom per element, as seen on the right in Figure 4.

6.2 Comparing space-time meshing strategies

In Section 4, we have seen how to advance the solution element wise on a tent pitched mesh. We now compare this approach to solving the full system on a Cartesian (in time) space-time slab. To solve the full system we use a block Jacobi solver. When comparing the timing of the two methods, we consider 4 different cases for the tent pitching approach, first solving the tents sequentially, and then solving them in parallel on 6, 12, and 24 threads. For this comparison, we choose a quasi-uniform mesh of the unit square in space and the final time equal to the mesh size, i.e. one CFL-conforming time step on the Cartesian mesh. For the $p$-version comparison in Figure 5 on the top we fix $h=0.04$, and in turn, for the $h$-version comparison on the bottom we fix $p=3$.

The results, in Figure 5 on the left, show that the error between the two mesh types differs only slightly. On the right in Figure 5, we compare the runtime of solving the full system on a Cartesian mesh, with solving the tents sequentially (on 1 thread), and solving them in parallel. Sequential tent pitching is about one magnitude faster. Moreover, we can investigate the effects of parallelising the computations. Using more threads only gives an advantage for small enough mesh sizes, as we are only able to solve independent tents in parallel.

6.3 Choice of spatial basis functions

As we have seen in Section 3, the recursion formula (3.1) for the derivation of the Trefftz basis functions, can be initialized with an arbitrary choice of polynomial basis functions in space variable only.

In the following, we compare three different choices for the initial polynomial basis functions: monomials, Legendre, and Chebychev polynomials. In all cases, the basis functions are shifted to the center of each element and scaled, as described in Remark 3.1. We compare them in 1+1 dimensions, on the space-time unit square. The mesh considered is the tent pitched mesh shown in Figure 3. The problem is solved globally using formulation (LABEL:eq:TDG). The results in Figure 6 show that all choices behave the same for low degrees. However, for higher degrees, Legendre and Chebychev polynomials fail to approximate the solution, due to the bad conditioning of the system matrix, compared to the monomials. The good properties of the two sets of basis functions do not carry over when developed in the recursion.

6.4 Tent pitching in 2 and 3 space dimensions

As discussed in Section 4, we solve elementwise, and in parallel. For this example, we choose as a spatial domain $\Omega$ the unit square and the unit cube. The initial quasi-uniform spatial mesh consists of triangles or tetrahedrons of maximal size $h$. We then use tent pitching in $2+1$ and $3+1$ dimensions, until the algorithm stops at time $T=1$, where we compute the error. The results of this are shown in Figure 7.

In Figure 7 on the left, we plot the error in terms of $h$ for different values of polynomial degree $p$ of $\mathbb{W}^{p}(Q)$. As typical for DG methods, in the case of a regular enough solution, we observe superconvergence, with the rate $\mathcal{O}(h^{p+1})$, outperforming the expected $\mathcal{O}(h^{p+1/2})$, see [20, Thm. 6.19]. We also consider convergence in terms of degree $p$ of the Trefftz space $\mathbb{W}_{p}(\mathcal{T}_{h})$, and report the results in Figure 7, right plots. For our analytic solution, we can observe exponential convergence.

6.5 Dissipation of energy

For smooth enough functions $(w,\bm{\tau})$ the energy at a fixed time $\hat{t}$ is given by

[TABLE]

In [20], the method (LABEL:eq:TDG) was shown to be dissipative, which we can also observe in numerical examples. We test on a model problem with analytical solution

[TABLE]

on the domain $[0,1]\times[0,T]$. We solve using the tent pitching algorithm.

The space mesh considered is a uniform partition of the interval $[0,1]$ into 5 elements. We measure the relative error in the energy given by

[TABLE]

In Figure 8 on the left, we can see that the error in energy increases in time. As the energy of the analytical solution is constant, and we are plotting the error without absolute value, we deduce that the energy of the numerical solution is decreasing in time. The results actually suggest that the energy of the numerical solution decreases linearly in time. In Figure 8 on the right, we compare the error in the energy at three different times $T=100,1000,10000$, plotting it against the degree of Trefftz polynomials in a range from 2 to 8. We observe exponential convergence for increasing order. Furthermore, greater times $T$ seem to affect the error only by a multiplicative factor.

6.6 Non-uniformly refined spatial meshes

Now that we have verified the convergence of the Trefftz-DG method with tent pitching initialized on quasi-uniform spatial meshes, we test the advantages of the method on a non-uniformly refined spatial mesh.

In this test, the refinement is applied to resolve a singular solution at the reentrant corner of an L-shaped domain, given by $\Omega=[-1,1]^{2}\setminus([0,1]\times[-1,0])$. The mesh refinement strategy used takes the diameter of a spatial mesh elements $K$ as

[TABLE]

where $r$ is the distance of $K$ to the reentrant corner, fixing a minimal mesh size of $h_{\min}=h_{\max}^{1/\mu}$. Motivated by the theoretical results in [1], we choose $\mu=\frac{1}{3}$.

We consider a model problem with solution given, in polar coordinates, by

[TABLE]

where $J_{\nu}$ denotes the Bessel function of the first kind. We consider $\nu=2/3$, so that $\nabla U$ is singular at the origin. We solve up to time $T=1$ for $a=10$. To avoid numerically integrate the singularity, we use the method to reconstruct the second order solution $U_{hp}\in\mathbb{U}^{p}(Q)$, introduced in Section 5, to measure the error given by $\|U(\cdot,T)-U_{hp}(\cdot,T)\|_{L^{2}(\Omega)}$.

The comparison between results obtained with uniform and non-uniform mesh refinement are shown in Figure 9 for Trefftz functions of degree $p=3$. We compare the two different meshing strategies by plotting them against (global dof)-1/3. For the uniformly refined meshes, the convergence rate is bounded by the smoothness of the solution $U\in H^{5/3-\varepsilon}(Q)$, for $\varepsilon>0$. We observe a convergence rate of $\mathcal{O}(h^{1})$. Using the non-uniformly refined meshes, we are able to recover optimal convergence for the third order Trefftz polynomials, as seen in Figure 9.

Table 1, gives a closer look on some of the properties already visualized in Figure 9 and also shows the runtime (in seconds). For the computations we used 24 threads. In Figure 9 on the bottom right, we compare the run time with the degrees of freedom. We observe that the uniform and the non-uniform mesh take about the same time for comparable numbers of degrees of freedom. Thus, no significant locking, due to the spatial refinement, occurs.

6.7 Wave propagation in an heterogeneous material

In the following example we investigate the reflection of a wave at an interface of two different materials. This experimental setup was also performed in [16, 2]. We consider the space-time domain $Q=[0,2]^{2}\times(0,1]$, and problem (LABEL:eq:TDG) with homogeneous Dirichlet boundary conditions. The wavespeed is the piecewise constant function given by

[TABLE]

As initial condition, we take a Gaussian wave given by

[TABLE]

where we choose $\bm{x}_{0}=(1,1)$ and $\delta=0.01$. The computations are performed with polynomial degree $p=4$.

Snapshots of the solution are shown in Figure 10. In the Snapshots, the right part of the domain has spatial mesh sizes up to $0.03$, whereas in the left part we choose as spatial mesh size of $0.01$, in order to better capture the steeper wavefront in the slower traveling material. First, we see that the initial condition unfolds in the left homogeneous part of the medium. At $T=0.2$, the wave crosses over into the material with higher wave velocity. In the next snapshot we can see that the wave splits into a part traveling to the right with a higher velocity and shallow wavefront, and a part reflected at the interface traveling backwards to the left. Finally, at the time $T=0.4$, we can also observe the weaker Huygens wave, which traveled parallel to the interface, before traveling back towards the left.

In Figure 11 on the left we present a sketch of the actions described above, also indicating a region where we measured the output

[TABLE]

The domain of measurement was chosen $\Omega_{C}=[1-\varepsilon_{C},1+\varepsilon_{C}]\times[0.25-\varepsilon_{C},0.25+\varepsilon_{C}]$, with $\varepsilon_{C}=2^{-7}$. The measurement over time is presented in Figure 11 on the right and shows that we are able to distinguish the three incoming waves. We can see the very weak Huygens wave arriving first, followed by the initial wave and the reflected one.

7 Conclusion

We have presented implementational aspects and numerical results for the Trefftz-DG method for the acoustic wave equation, originally presented in [20]. The implementation in NGSolve was used to solve on Cartesian (in time) meshes and tent pitched meshes in up to 3+1 dimensions, with varying mesh sizes, polynomial degrees, and wavenumber.

The $h$-convergence rates were shown to comply with the analytic results, showing superconvergence for analytic solutions and limited rates in the case of solutions with insufficient regularity. In the latter case, we were able to recover optimal convergence using non-uniform meshes. For analytic solutions we observed exponential convergence rates in the polynomial degree $p$.

Possible developments include the extension to the case of electromagnetic waves (Maxwell’s equations).

Acknowledgment**.**

This work has been supported by the Austrian Science Fund, grants no. W1245 and F65.