High-order space-time finite element schemes for acoustic and viscodynamic wave equations with temporal decoupling

We revisit a method originally introduced by Werder et al. (in Comput. Methods Appl. Mech. Engrg., 190:6685–6708, 2001) for temporally discontinuous Galerkin FEMs applied to a parabolic partial differential equation. In that approach, block systems arise because of the coupling of the spatial systems through inner products of the temporal basis functions. If the spatial finite element space is of dimension D and polynomials of degree r are used in time, the block system has dimension (r + 1)D and is usually regarded as being too large when r > 1. Werder et al. found that the space-time coupling matrices are diagonalizable over for r ⩽100, and this means that the time-coupled computations within a time step can actually be decoupled. By using either continuous Galerkin or spectral element methods in space, we apply this DG-in-time methodology, for the first time, to second-order wave equations including elastodynamics with and without Kelvin–Voigt and Maxwell–Zener viscoelasticity. An example set of numerical results is given to demonstrate the favourable effect on error and computational work of the moderately high-order (up to degree 7) temporal and spatio-temporal approximations, and we also touch on an application of this method to an ambitious problem related to the diagnosis of coronary artery disease. Copyright © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.


INTRODUCTION AND MOTIVATION
In 2001, Werder, Gerdes, Schötzau and Schwab, [1], formulated a space-time FEM for the heat equation that employed a discontinuous Galerkin (DG) method in time. This was not the first such formulation (see e.g. [2]), but one particularly interesting feature of their work demonstrates that the temporal system can be diagonalised, thus making reasonably high-order finite element time discretizations feasible. Also, there are recent results [3,4] (see also the much earlier [5]) showing that high-order spatial discretizations for wave equations are highly desirable in terms of the control of dispersion errors, and it is against this background that we aim here to extend the work of Werder et al. to facilitate high-order space-time discretizations for second-order hyperbolic wave equations. Our goal is to present a framework for useful and practical high-order space-time FEMs for a collection of common and important linear wave equation problems. We remark also that while high-order time stepping has been around for some time, by using Runge-Kutta methods for example, the approach in the succeeding text allows for the entire space-time discretization to be placed in a variational setting. This allows the arsenal of very powerful functional analytic tools relating to stability and error estimation to be deployed for the numerical simulation of wave propagation.
We deal firstly with an abstract formulation of the wave equation in Section 2 and describe the time discretization and its diagonalisation. We then apply this abstract framework to the specific examples of the acoustic wave equation and the equations of elastodynamics and viscodynamics in Section 3. The wave equation is dealt with rather swiftly in Subsection 3.1, and we move on to linear elasticity and Maxwell-Zener (hereditary) viscoelasticity in Subsection 3.2 and, finally, the decoupling procedure in the presence of both long-memory and short-memory viscoelasticity is described in Subsection 3.3.
A collection of numerical results is given and discussed in Section 4, where we aim to demonstrate the convergence rates that are achieved in the natural norms associated with the wave equation. Our presentation is entirely practical in that we do not aim to establish theorems on error bounds but rather to demonstrate the practicality of high-order finite element time stepping. Our motivation for adopting this methodology is to simulate waves in biotissue with a specific application related to the diagnosis of coronary artery disease. This application and the challenges it presents are outlined in Section 5. Section 6 then concludes with some observations. We would like to point out that this is a shorter version of the long report in [6]. That report contains a much more extensive set of numerical tests.
The approach to discretization taken in the succeeding text utilises a normalized Legendre basis to effect the DGFEM in time discretization while in space we report here on the standard continuous Galerkin (CG) spatial discretization. The longer report, [6], contains results also for the case where the continuous spectral element method (SEM) using Gauss-Lobatto integration and nodes is used for space discretization.
Irrespective of the use or not of the spectral FEM, we shall see in the succeeding text (in (8) and (9)) that the decoupling procedure produces a set of boundary value problems that need to be solved for each time interval. Any suitable method could be employed here, such as the standard Galerkin procedure, a 'blended procedure' as in [7], or indeed any of the many other tools that have been developed to 'solve' elliptic problems (such as, if applicable, exact solutions and asymptotics).
To provide some more context for what follows, we recall that the idea of using Galerkin finite elements in space dates back to [8] and [9] and, of course, is now a method of choice in the vast majority of cases involving elliptic operators. However, although given much less exposure, the idea of using Galerkin discretizations in time is also not new (e.g. [10][11][12]) but it seems clear that these methods never really caught the imagination of users and producers of codes until much more recently. In fact, arguably, it was the work of Eriksson and Johnson, see [13], on adaptive spacetime formulations for parabolic problems that seemed to have revived interest in these methods. The formulations in [14][15][16][17][18] are closely related but here we stay much closer to the crisper formulation given by Johnson in [19]. We consider only DGFEM in time in this effort specifically because we are building on the work in [1], although we note that there also exist CG methods for time discretization (e.g. [20,21] and [22]). These are important in their own right, due mainly to their stability and energy conservation properties, and our early results for the heat equation in [23] suggest that a 'CGFEM-in-time' could be developed for the wave equation. We leave this for another time.
We close this introduction by recalling one very basic reason why high-order schemes are useful. Suppose we compute up to time T using N 1 1 time steps. If piecewise, polynomial degree r 1 is used for the approximation, then on each step we can expect the error to be of the order N d 1 133 computational work needed is then N 1 d 1 =N 2 d 2 D N 1 1 1 as N 1 and/or become large and so we conclude that, with sufficient solution regularity, higher-order schemes are capable of providing higher-fidelity solutions than lower-order schemes for the same amount of computational work. We attempt to illustrate this later for the 2D examples in Section 5.

AN ABSTRACT FRAMEWORK FOR DECOUPLED DISCONTINUOUS GALERKIN FEM IN TIME
Our notation is standard and, apart from the preliminaries that follow, is introduced where necessary. Let the spatial domain of interest, , be a time-independent open-bounded polytope in R d for d = 1, 2 or 3. We assume that the boundary, @ , is partitioned into ¹ D ; N º (also time independent) and assume that Dirichlet boundary values are given on the closed set D with Neumann boundary values specified on the open (and possibly empty) set N . As usual, we will insist that D \ N D ¿ and D [ N D @ . We note that, in general, there exist problems for which D may have surface measure zero, but here we insist that meas @ . D / > 0. Also, the unit outward normal vector to N will be written as O n. To deal with the time dependence, we set I WD .0; T and will usually use overdots to denote partial time differentiation.
Before getting to the specific formulations of the later sections, we firstly introduce an abstract formulation along with a semidiscrete time discretization using the DGFEM. To this end, we let V ,! H ,! V 0 be a Gelfand triple of reflexive Hilbert spaces (dense and continuous embedding, see, e.g. Wloka [24]) and denote the inner product and norm on H by . ; /W H H ! R and k k D . ; / 1=2 . We deal in this section with the standard abstract form of second-order wave equation problems in order to fix ideas and notation. Later, in Section 3, we recall some concrete applications and also some damping terms. These are important for the soft-tissue applications that we have in mind (see later in Section 5). Our goal here is to formulate a high-order finite element time discretization that is rendered practical by decoupling. However, we point out once again that such a scheme will only result in higher-solution quality if there is sufficient solution regularity.
Let aW V V ! R be a symmetric and V -coercive bilinear form and assume that, for almost every t 2 I , we are given data L.t /W I ! V 0 . Introducing w WD P u and suppressing spatial dependence for clarity, we consider the abstract problem of finding a smooth map uW I ! V such that, where h ; iW V 0 V ! R denotes the duality pairing arising from the continuous extension of . ; / H to V 0 V and M u; M w are initial data. For convenience and the avoidance of unimportant constants, we henceforth consider the space V as equipped with the scalar product a. ; / and the induced norm k k V D a. ; / 1=2 . In this, % is a known positive constant that in solid dynamics represents the material's mass density. Also, and as is usual for time-dependent problems, for a Banach space .B; k k B /, we define the L p .0; t I B/ norms by kvk L p .0;t IB/ WD k kv. /k B k L p .0;t / . The DG-in-time formulation that follows is identical to that given by Johnson in [19]. To give a little background, we recall that in [19] the time dependence of the unknown displacement and velocity is approximated by piecewise linear functions that are allowed to be discontinuous at the time nodes. Because each function has a first-time derivative taken, we then expect to see 'delta functions' at these nodes that, when integrated (as in a Galerkin formulation), produce jumps in function values. That is the reason for the jump notation that is defined and used in the succeeding text. A full and rigorous derivation of the weak formulation that follows is beyond our scope here but can be accomplished by appealing to the definition of a weak derivative or to the Theory of Distributions (e.g. [25]).
Discretizing the time interval so that 0 D t 0 < t 1 < < t N D T , defining the time step k n WD t n t n 1 and setting I n WD .t n 1 ; t n / we define, as usual, the jump notation, The semidiscrete DGFEM approximation of (1) is then: for each n D 1; 2; : : : ; N in turn, find .U; W /j I n 2 P r .I n I V / P r .I n I V / such that Z with the understanding that U 0 WD M u and W 0 WD M w. Here, for each n, we use P r .I n I X/ to denote the space of polynomials of degree r on the time interval I n with coefficients in the target space X . The target space is omitted when X D R, and we note that r could be n-dependent. Choosing # D W and D U gives, and because we arrive at Taking L D 0 and summing over n D 1; ; N then gives the basic stability estimate, and reminds us that this scheme is dissipative. This, of course, is a general statement and holds for all polynomial degrees. It is also useful at this stage to recall Claes Johnson's result in [19] which focussed on Galerkin schemes by using discontinuous piecewise linears in time and continuous piecewise linears in space for the standard problem R u r 2 u D f subject to homogeneous Dirichlet boundary data. Although the presentation of the a priori error bound, there is rather terse it seems from the comments in the succeeding text and in the introduction that we can expect the following, 135 generalization of these error bounds to higher-degree polynomials (although we might expect this to be straightforward). We now move on to the specifics of the implementation. Let ¹ i W i D 0; 1; ; rº be a basis for P r .I n / and introduce the ansatz forms of the approximations to u and w on I n as, where ¹U 0 ; U 1 ; : : :º; ¹W 0 ; W 1 ; : : :º Â V . Replacing each of #.t/ and .t/ with i .t /# for i 2 P r .I n / and # 2 V in (3), we obtain firstly, where each holds for all # 2 V and for each i 2 ¹0; 1; : : : ; rº. Define matrices via, where i indexes the rows, and then by choosing our basis functions as the image under the linear map from OE 1; 1 to I n of the normalized Legendre polynomials, we conclude easily that 2M D k n I. It seems useful at this point to recall the remark in [1] that these matrices are hierarchical, and also that they can be precomputed just once on a reference length and then reused for subsequent computations. However, the main point and motivation for us here is that Werder et al. in [1] report that A is diagonalizable over C for all polynomial degrees of practical interest. This in particular means that the computational linear algebra associated with these objects is relatively cheap, and it also allows us to write D D Q 1 AQ D p 1 rC1 y where p y indicates a diagonal matrix of pairwise complex conjugate eigenvalues and where Q has complex entries. So, proceeding step by step and using the abbreviation † j to mean the summation P r j D0 , our where are known from data and the previous time step. Defining ¹Y q º and ¹Z q º as the unique solutions of and taking linear combinations using the rows of R WD Q 1 then yields, The algorithm is now clear. At each time step, we solve r C 1 boundary value problems for the complex unknowns ¹Y i º r i D0 and then update with We recall that for the SEM, the mass matrices arising from the spatial discretization of these systems are diagonal, and note also that because 2 i is complex we will either need to introduce complex arithmetic in the solver or treat the complex system as a two-by-two-block real system.

Remark 2.1 (non-zero Dirichlet data)
Although not included in the formulation in the previous text, some of the examples in the succeeding text use non-zero Dirichlet boundary data. These data are imposed by time-projection on to the temporal basis at each desired x 2 D to determine pointwise boundary values for each of the decoupled systems. Specifically, suppose that we want to set the known Dirichlet value w D .t / Á P u D .t / during the time interval I n . We assume the ansatz w D .t / D P r j D0 j .t /w D j (with the superscript D denoting 'Dirichlet') and then determine the w D j by the projection: We then do the same for the U D j values at the boundary and we modify (6) (10) where W n 1 D W n 1 I C W n 1 D is now used in the definition of F i , and the superscript I denotes 'interior'. The W D j and U D j on the right are then spatially interpolated along the boundary by augmenting the spatial basis to incorporate the Dirichlet boundary nodes (with zero values in the interior), and then 'mass' and 'stiffness' (rectangular) matrix-vector contributions to the load are formed from these new inner products on the right hand side of (10). Because we derive the boundary values for w from P u, they cancel each other out in (7) and so no change is needed there. The remainder of the implementation is as described earlier but with the 'interior' functions being the unknowns.

Remark 2.2
We could also define the semidiscrete DGFEM approximation of (1) as: for each n D 1; 2; : : : ; N in turn, find .U; W /j I n 2 P r .I n I V / P r .I n I V / such that Z 137 The advantage of this is that we could deal with pure Neumann problems wherein a. ; / may no longer be coercive. Numerical tests show that this scheme works well, but we are not able to prove any stability estimates for it. For that reason, we do not consider it further here, but note that the time decoupling method described in the previous text could be applied to this formulation also.

SPECIFIC APPLICATIONS
In this section, we apply the foregoing to the specific examples of the wave equation and elastodynamics. We also enhance the formulation to take viscoelastic damping effects into account.

The acoustic wave equation
In this subsection, we consider the specific example problem where we seek u such that, Defining a.w; v/ WD c 2 rw; rv , where . ; / is the L 2 . / inner product, and setting V WD ¹v 2 H 1 . /W vj D D 0º the weak form of the problem is: find uW I ! V such that, where hL.t /; i WD .±; / C .g; / N . Notice that we overload . ; / in a completely standard way where no notational distinction is made when the arguments are scalar-valued or vector-valued. The semidiscrete finite element formulation in (3) applies without change and so on the n-th time step, of width k n D t n t n 1 say, and for i D 0; : : : ; r, we have to solve the boundary value problems given by: where Ä i WD 1 4 .k n = i / 2 , and perform the update The spatial discretization has already been described in general terms earlier in Section 1, and so all that remains is to demonstrate the behaviour of the scheme. This is carried out later in Subsection 4.1 for both 1D and 2D test problems. In those experiments, we will have either c 2 D E 0 =% or c 2 D G=% where E 0 and G represent stiffness moduli and % a mass density (Section 5).

Linear elastodynamics and viscodynamics
In this section, we suppose that represents the interior of a linear viscoelastic compressible body. This body is acted upon by a system of body forces ± WD .± i .x; t// d i D1 for x WD .x i / d i D1 2 and t 2 I and on the open (and possibly empty) set N there is prescribed a system of surface tractions g WD .g i .x; t// d i D1 for x 2 N and t 2 I . The displacement from equilibrium resulting from the action of the applied forces ± and g is denoted by u D .u i / d i D1 W I ! R d and, in this linear theory, the deformation is described by the strain tensor " WD ." ij / d i;j D1 given by, We assume also that t D 0 is a reference time such that u D 0 for all t < 0. Newton's second law of motion with boundary conditions gives for each i 2 ¹1; : : : ; d º that, with summation implied and where WD . ij / d i;j D1 is the symmetric stress tensor. We close this problem by recalling the standard literature on viscoelasticity (e.g. [26,27]) and introducing the following linear hereditary viscoelastic constitutive relationship between stress and strain, where, as usual, we omit the x dependence and, as in the previous text, sum over repeated indices. In this, C and D.t / are fourth-order tensors with C related to the Kelvin-Voigt model of viscoelasticity and D.t / stemming from the Zener and Maxwell models. The Hooke's law tensor D.t / WD D ij kl .t / d i;j;k;lD1 is a fourth-order stress relaxation tensor satisfying the following symmetries: However, we do have D ij kl .t / D D klij .t / for t D 0 and t D 1 in general, and for all t for isotropic materials (see, e.g. [28, equations (1.10), (2.62)]), and the components of D can be assumed to be (a.e. in ) of class W 1 1 in t although we will not need this level of generality. In addition, because D.0/ measures instantaneous linear elastic response, we follow Hooke's law and assume positive-definiteness: ij kl D ij kl .0/ > 0 a.e. in for all non-zero symmetric second-order tensors .
In what follows, we will assume a much simpler version of this constitutive law. Firstly, we assume that the material is synchronous so that D.t / is replaced by '.t/D, where D is now constant in time, and ' is a stress relaxation function. Following classical theory, we assume the Prony series form, where ' q > 0 for q 2 ¹0; 1; : : : ; N ' º, q > 0 for q 2 ¹1; : : : ; N ' º, and we normalize so that The case, ' 0 D 0, corresponds to a viscoelastic fluid in the sense described by Golden and Graham in [26], whereas ' 0 > 0 gives a solid. The second assumption is that the material is homogenous and so, in particular, is isotropic. This means that D can be described by just two independent Lamé coefficients denoted, usually, by D E=..1 C /.1 2 // and D 2G D E=.1 C /, where E is Young's modulus, G is the shear modulus and is Poisson's ratio. Theoretically, 2 . 1; 1=2/ but because we are dealing only with compressible non-auxetic materials, we have 2 .0; 1=2/ and both and are well defined. The action of D is now given specifically by D ij kl " kl .u/ D r uı ij C " ij .u/, and the constitutive law can be written in the form, where for q D 1; 2; : : : ; N ' , are internal variables (see, e.g. [29][30][31][32]) and satisfy, To give a weak formulation of this problem, we firstly define the product Hilbert spaces, H s . / WD H s . / n , for s D 0; 1; 2; : : :, with inner products given for all , v 2 H s . / by These spaces have the natural norms k k s WD p . ; / s and, of course, L 2 . / Á H 0 . /. We also use the (symmetric second-order) tensor-valued L 2 space, and, noting the essential boundary condition, we define the test space as, After integration by parts (see e.g. [33] for details), we then arrive at the weak problem of seeking a smooth map, uW I ! V , such that, This completes the statement of the basic problem. The next section discusses the extension of the decoupling procedure to incorporate the viscoelastic terms that have just been introduced.

Decoupling in the presence of viscoelastic damping
Now that we have defined the viscoelastic models, we introduce them into the DG formulation as given earlier by (3) and obtain the following problem. For each n D 1; 2; : : : ; N in turn, find U; W; U 1 ; U 2 ; : : : ; U N ' j I n 2 P r .I n I V / P r .I n I V / P r .I n I V / N ' such that, 8# 2 P r .I n I V /; 8¹Á q º N ' qD1 Â P r .I n I V / and 8 2 P r .I n I V /: The decoupling is now a straightforward extension of the method described earlier and by similar means, and with the same notation for the temporal basis, we obtain firstly, and thirdly, i;n 1 a .U q / n 1 ; Á q Ã D 0: In our earlier notation, these are simply and with easily derivable right hand side functionals F, G, G 1 ; : : : ; G N ' . The decoupling is now carried out exactly as before.

NUMERICAL EXPERIMENTS
In this section, we give a selection of results from some computational experiments. Our main goal is to illustrate the estimated convergence rates, for the standard norms, which are achieved by this scheme. The examples have been chosen in order to try and cover all of the main possibilities. For example, we can have mixed or Dirichlet boundary data along with, in each case, various combinations of elastic and viscoelastic effects. In Subsection 4.1, we focus on the scalar-wave equation as described earlier in Subsection 3.1. Then in Subsection 4.2, we give some similar results the viscoelasticity problem with the scheme described earlier in Subsection 3.3. As for practical matters, we remark firstly that the enforcement of initial data was, in all cases, carried out by interpolation rather than projection and, secondly, that all integrals were computed by using Gauss-Legendre 141 rules of high-enough order to render the contribution of the 'variational crime' negligible as compared with the errors that we are interested in. To detail this recall that the Gauss-Legendre n point rule is exact for polynomial integrands of degree 2n 1, and we ask that this quadrature rule be exact for three orders higher than the maximum degree that occurs. This means that for an approximation based on polynomials of degree r, we require 2n 1 D 2r C 3 and so n D r C 2. This is the order of rule used (in space and time -whenever necessary) in our 2D results. As alluded to in the introduction, we have also computed results for spatial discretization using the SEM where Gauss-Lobatto rules are used to build the system matrices with nodes placed at (the tensor product of) those same Gauss-Lobatto nodes. The Galerkin method for which results are given in the succeeding text used the same node locations. This was purely for ease of implementation where the same code could perform both the Galerkin and the spectral computation.

The acoustic wave equation
We consider several examples, focussing first on a one-space dimensional problem in order to best reveal the (computed) temporal convergence rates. Specifically, for examples 1, 2, 3 and 4 that now follow, we consider a problem that involves different combinations of the Kelvin-Voigt and Maxwell-Zener damping. The problem, with x-dependence suppressed, is where D .0; 1/ with T D and all data chosen so that u D x.x 1/ cos.t / is the exact solution. Homogeneous Dirichlet boundary data are imposed at each end of , and we use both the Galerkin and spectral spatial discretization. In these computations, we used 12 point Gauss-Legendre integration to compute the spatial inner products 'exactly' and, for the Galerkin spatial approximation, we used a two equal-element mesh for and piecewise quadratics.
For the DG time discretization, we use piecewise polynomials of degree 1 through 7 and also compute using a Crank-Nicolson (CN), or trapezoidal, discretization of (19) and (20) to give a comparison with a standard and well-known scheme. In that, 'CN method', we put w D P u and then use central differences across a time step for the time derivatives and time nodal-averages for the remaining terms (this is a very standard discretization, see [34] for a similar approach).
We show only a small selection of the computed results here. The full set are available in the extended report [6].

Example 1 (time error; undamped)
In (23), we choose the density as % D 1010 kg=m 3 , the stiffness modulus as E 0 D 58 kPa, damping modulus as E 1 D 0 Pa s and '.t/ D 0. To save space, the results are not shown because they are very similar to those in Example 4 in the succeeding text.

Example 2 (time error; Kelvin-Voigt)
This is exactly as mentioned previously for Example 1 but with E 1 D 30 kPa s. The results for the computed energy error bound given by (4) are shown in Figure 1 for the Galerkin error. The SEM gave a similar picture.

Example 3 (time error; Kelvin-Voigt and Maxwell-Zener)
Similar to Example 2 but with '.t/ D 1 C e t =0:05 =2. The results are not shown because they are very similar to those in Example 2 in the previous text.

Example 4 (time error; Maxwell-Zener)
Again, similar to Example 2 but with E 1 D 0, and '.t/ D 1 C e t =0:05 =2. The results for the computed energy error bound given by (4) are shown in Figure 2 for the Galerkin error. The SEM gave a similar picture.   triangles'. These indicate the slopes that correspond to energy error bounds of size O.k p / for p 2 ¹2; 3; 4; 5; 6º, for the fan on the right and p 2 ¹7; 8; 9º, for the fan on the left. For the work-error plots, we define work as the product of the number of time steps with the number of matrix solves required per time step. For the decoupled DG approximation using polynomials of degree r, this latter quantity is of course r C 1. For comparison purposes, we also show the performance of the CN method, recalling that it requires only one real-arithmetic solve per time step. Assuming that real-arithmetic multiplications are four times faster than complex, that realarithmetic additions are twice as fast as complex and that banded-matrix (of dimension, D, say) inversion is of quadratic complexity with an equal mix of products and sums, we can estimate the operation count of a real inversion as .D=2/ 2 C .D=4/ 2 D 5D 2 =16 and that of a complex inversion as 2D 2 . The complex inversion is therefore around 32=5 D 6:4 times more expensive. In producing the work-error data for the figures, we took one complex inversion as the basic work unit (per time step) and applied the scaling of 5=32 to the CN figures to make the comparison fairer. Actual raw wall-clock timings from the Example 1 run produced a real scalings of around 1=6 to 1=6:5, so our approach appears to be well-founded. Of course, for larger and more realistic problems, an iterative solver will almost certainly be required and the work comparison becomes a much more delicate issue.
The results for the 1D problems in Examples 1, 2, 3 and 4 show the behaviour of 'time error' in the computed solution for four cases covering elasticity and mixtures of viscoelastic formulations. Copyright

143
Because Examples 1 and 4 were similar, and Examples 2 and 3 were also similar, we can infer that the presence of Maxwell-Zener viscoelasticity does not alter the convergence properties of the scheme away from that displayed by the standard wave equation. The addition of a Kelvin-Voigt term, however, does seem to negatively affect the convergence behaviour but without a complete a priori error analysis to hand it is not clear why or how.
In particular, for Examples 1 and 4, from Figure 2, we can infer CN convergence with energy error like k 2 , DG(1) like k 3 , DG(2) like k 5 and DG(3) like k 7 . It seems safe to conjecture a DG(r) energy error of order O.k 2rC1 /.
On the other hand, for Examples 2 and 3, from Figure 1, we can infer that the CN scheme converges with energy error like k 2 and the linear DG scheme with energy error like k 3 . These are as expected. However, the quadratic and cubic DG schemes seem to converge like k 4:5 and k 6 , whereas the quartic and quintic seem both to be between k 7 and k 8 . The sextic and septic schemes appear to converge with at least k 9 . The reason for these unexpected behaviours is not clear although we suspect that the precision of the computations is inadequate for these high-order approximations. It does, however, seem safe to conclude that convergence is not so rapid when a Kelvin-Voigt term is present.
The accompanying work-error plots seek to illustrate the efficiency associated with higher-order methods. However, they need to be interpreted with care because it would be unlikely to have only 'time error' in a problem of real interest. They are included merely to remind us that the extra linear algebra workload must be borne in mind when decoupling in time. These results were for the Galerkin-in-space method. The results for the SEM were similar in both cases.
We turn now to problems in two space dimensions and an illustration of space-time error behaviour. The physical data are based on the 'real life' problem that we discuss later in Section 5, and the spatial domain is a scaled-up version of the one that we later describe.
For this 2D calculation, we chose the domain D ¹0:005 < x < 0:15 and 0 < y < 0:3º with T D 2 and the problem % R u Gr 2 u D ± with % D 1010 kg=m 3 and G D 58 kPa. We took N uniform time steps and, to demonstrate the computed space-time convergence rate, we meshed the domain with an N by N grid. These numerical results use DG-in-time with piecewise polynomials in space and time of degree up to 7, with the same degree being used in both space and time except when indicated. The CN computations are based on piecewise linears for the spatial approximations, with quadratics being used in the bottom row of the figures-as explained more fully in the first example. Our motivation for mixing these degrees was to check whether there was any numerical evidence of a different convergence rate 'in space' than 'in time'. The bottom rows show that in fact there is -we can expect at least one degree higher in time for the L 2 . / error in both u and ru. In the convergence graphs that follow, the triangle fans indicate convergence rates from zero to nine in half-steps.

Example 5 (space-time error for the wave equation, Dirichlet BC's)
Here, Dirichlet boundary data are assumed on the whole of @ and, ±, the boundary and the initial data are chosen so that the exact solution is u D cos.2 t / sin.40 x/ cos.30 y/. The results are plotted on the left of the top row in Figures 3, 4 and 5 for the Galerkin-in-space scheme for errors in u.T /, ru.T / and w.T /, with the estimated convergence rates tabulated on the right. (Notice that a zero value of L 2 error for w is reported for the CN scheme when N D 2 -this is no more than a numerical anomaly resulting from exact boundary data and a very coarse mesh.) These are based on the computed errors for N (first column) and for 2N (the results for N D 4 use those for N D 2, the row for which is not shown). The middle row of the figure shows the result of computing with one degree higher polynomials in time than in space, whereas the bottom row shows the result of using one degree higher in space than in time. The CN results are of course unaffected in the middle row. Analogous results for the spectral element scheme in the case where equal degrees are used have also been obtained and are contained in [6].  With Examples 5 and 6 for the scalar wave equation, we can now discuss space-time convergence. Firstly, for Example 5, we infer from Figure 3 for the Galerkin-in-space scheme, with the last term being suggested by the bottom table in the figure.
In a similar way, Figure 4 suggests that ru.t N / rU N L 2 . / D O h d C k d C1 , although the last term is not so clearly apparent in this case. The results for w D P u are not so useful. The curves in Figure 5 quite clearly show rapid convergence as the order increases but the tables do not give useful figures. The convergence behaviour seems too erratic to be able to draw useful conclusions. Our longer report, [6], contains the analogous results for the SEM and here we simply remark that the same conclusions can be drawn for that as for the Galerkin scheme. Example 6 has different types of boundary data. Figures 6, 7 and 8 for the Galerkin-in-space scheme and the results for the spectral element scheme in [6] suggest essentially the same behaviour as for Example 5.

Elastodynamics and viscodynamics
In this subsection, we give a numerical demonstration of this scheme for the elasticity and viscoelasticity problems described earlier in Subsection 3.2. Also, unless stated otherwise, the set-up for the numerical results is exactly as for the wave equation as described in the previous section.

Example 7 (space-time error for viscodynamics, Mixed BC's)
For this 2D calculation, we choose the same domain as previously mentioned, D ¹0:005 < x < 0:15 and 0 < y < 0:3º but now with T D 0:5 and the problem given by (19) and (20). We set % D 1010 kg=m 3 , E D 167 kPa, Poisson's ratio D 0:44 and include viscoelastic effects by choosing ' 0 D 0:2, ' 1 D 0:8 and 1 D 0:05 in (16). The load, ±, boundary and initial data are chosen so that the exact solution is .u 1 ; u 2 / T D .sin.2 x/ sin.2 y/ cos.2 t C =4/; cos.2 x/ sin.2 y/ sin.2 t =4// T and Dirichlet data were prescribed on bottom edge of @ with tractions prescribed elsewhere. The results are given in Figures 9, 10 and 11 for the Galerkin-in-space scheme and, as before, analogous results for the SEM are in [6]. These results are for the case where equal polynomial degrees in space and time were used.  We observe that Figures 9, 10 and 11 for the Galerkin-in-space scheme are compatible with the expectations of error rates indicated by Examples 5 and 6. For that reason, we did not pursue a more extensive numerical study of convergence in this case (again, see [6] for the SEM).
Although these rates of convergence are, broadly, in line with what we might expect, we should sound a note of caution. A single-discretization parameter, N , has been used which means that the mesh width and time steps are, in essence, the same. The mixing of h and k suggested by Johnson's results ((4) and (5)) will therefore not be clearly revealed in our study. Also, curiously, we did not clearly observe the semi-integer rate predicted by Johnson's bounds, although it is possible that this is 'hidden' in the quite ragged results for w.

APPLICATION: SHEAR WAVES IN BIOTISSUE-MIMICKING GEL
In this section, we shift our focus away from exact errors and manufactured solutions to a much more practical experimental set-up. We consider an annular cylinder of tissue mimicking agarose gel of height 0:0514 m and with inner and outer radii 0:00175 m and 0:027 m. Here 'height' indicates that the cylinder is up-ended so that the radial plane is horizontal.
This experimental rig actually exists in our laboratories (see the right of Figure 12) and is the first step in an ambitious project investigating the possibility of diagnosing coronary artery disease through computational mathematics. Briefly, plaque build-up in a diseased coronary artery causes a stenosis (see the left of Figure 12) that induces a disturbance in the downstream blood flow. It  for Ń D 0:0257 and with a frequency, f , as specified in the succeeding text. Here, is a 'ramp-up' function the role of which is to ensure that high-frequency noise associated with sudden starts does not pollute the signal -this ramp allows for the full amplitude of the traction to develop slowly over 50 full cycles. This is intended to mimic the experimental set-up.
The material data are taken as % D 1010 kg=m 3 , E 0 D 167 kPa and D 0:44 (the Voigt term was set to zero: C D 0). High quality estimates for the viscoelastic parameters have now appeared in [39] but, for the sake of illustration here, we simply took ' 0 D 0:2, ' 1 D 0:8 and 1 D 2 with N ' D 1 in (16). The spatial discretization used in the succeeding text is a SEM on, mostly, a 13 (radial) by 25 (axial) uniform mesh of bicubics supported on tensor products of the Gauss-Lobatto nodes. (We also ran some problems on a 25 by 50 mesh, and will return to this point later.) The DG time discretization used piecewise cubic polynomials (we refer to this as 'DG(3)' in the succeeding text) and all except for the spectral quadratures (the mass and stiffness matrices) use high-order (degree 12 in space and time) Gauss-Legendre rules to calculate the (space and time) definite integrals.
The imposed traction, g, in (24) simulates the wall shear disturbance that may arise from a diseased (stenosed) artery, represented by the centre line of the cylinder and, in order to illustrate the Copyright Figure 4b in [35] ©RSNA 2009), an image of a right proximal coronary artery stenosis obtained by X-Ray angiography after injecting a contrast medium through a catheter passed into the aorta and thence into the root of the coronary vessel. On the right, an illustration of our experimental rig referred to in Section 5. The middle section of the scaffold contains the gel phantom wrapped in film to prevent dehydration. The rod containing a moulded-in bead is visible emerging from the top and the unit at the bottom is the controllable 'shaker' that imparts the bead's vertically linear vibration. This bead is mimicked in our computational set-up by the traction in (24). The displacement of the rod and bead assembly is measured by an optical device at the top of the rig, removed in this picture for clarity.     Selections from the computed signals are shown for f D 250 Hz and a 25 50 mesh in Figure 13, and for a 13 25 mesh in Figure 14. For f D 1250 Hz, selected output is shown in Figure 15, for the 13 25 mesh, whereas for the 25 50 mesh, we show output in Figure 16. Results for f 2 ¹250; 500; 750º Hz are given in [6].
To discuss these results, we begin with the computations for the least challenging frequency, 250 Hz, where we find output results for a selection of time step sizes and for the 25 50 spatial mesh in Figure 13. The corresponding results for the 13 25 mesh are shown in Figure 14.
In a crude 'eyeball judgement', we can see that the DG(3) scheme seems to have converged to a practically useful level of accuracy, and with both meshes, at 500 time steps, whereas CN requires around 12000 for both meshes. This equates to 2000 complex system solves for DG(3) as against 12000 real solves for CN.
We can use this observation to infer that it is reasonable to assume that the DG solution for the largest N DG is much more accurate than the CN solution at all considered values of N CN . This suggests a method for approximating the error in these computed signals. For each of the frequencies, we take the DG solution corresponding to the largest N DG in use, call these Q u and N DG;max , and assume that it is 'exact'. We then approximate the error in the CN and DG solutions on the set, S, of discrete times that are common to the 'exact' computation and the one being assessed. With these notions, we define the relative radial (subscript j D 1) and axial (subscript j D 2) errors by, The results for the resulting approximate relative error max¹E 1 ; E 2 º are shown in Figures 17 and 18 for the 24 50 mesh and the 13 25 mesh. We can see that they are effectively identical and, of course, we caution against paying too much attention to the estimated relative errors associated with the larger values of N DG . We also recorded the run times using MATLAB's ® tic and toc wallclock timer commands and although this allows the execution times to be plotted against number of time steps, we consider it more interesting to plot error as a function of execution time. Copyright   Here, for the 25 50 computation, we see from the left of Figure 17 that to obtain a relative error of about 10 3 we need around 500 DG(3) steps and 48000 CN steps. The right hand plot tells us that, for these values, CN took about 400 min: about 10 times that of DG(3).
Because Figure 17 and Figure 18 are essentially the same (except, of course, for execution time), we can assume that we have in effect computed an 'exact solution' that gives us confidence in our findings at 250 Hz.
In [6], we give 13 25 mesh results for 500 Hz, 750 Hz and 1000 Hz but include only the results for the highest frequency we computed for, 1250 Hz, here. In that case, the error and timing data are shown in Figure 19 for the 13 25 mesh and in Figure 20 for the 25 50 mesh.
It seems clear that on the 13 25 mesh the DG(3) results have converged by 3000 time steps and CN by 288 000; whereas for the 25 50 mesh, we would estimate these figures at 4500 and 288 000. However, the 'converged plots' are not the same -there are significant differences in the wave envelopes for times later than 0.2. This is a clear indication that a finer spatial mesh is needed in order to properly assess the situation at higher frequencies. This study is ongoing but will have to wait for a software development phase to be completed: the execution times required by MATLAB ® for repeating all these tests on a, say, 50 100 mesh, are prohibitive. Copyright

CONCLUSIONS AND DISCUSSION
The conclusions relating to the specific numerical performance of the scheme have been included in the previous two sections. Here, to close, we would like to make just two more general remarks. Firstly, despite the need to introduce complex arithmetic, it seems clear that DG (3) and, by extension, DG(p), can offer much greater efficiency for certain types of wave equation problem. Linearity, solution smoothness and (temporally) constant coefficients are the obvious desirable properties, and we note that in engineering dynamics these are often satisfied. The extension and study of the technique beyond these assumptions seems worthwhile, as does a full and detailed a priori error analysis building on that in [19]. On the basis of the results in the previous text, we could conjecture a temporal convergence rate at the nodes corresponding to a O.k 2qC1 / error term (as for the DG scheme for parabolic problems in [13]) when the Kelvin-Voigt term is not present. It also seems clear that the addition of Maxwell-Zener viscoelastic damping does not affect the convergence rates.
Secondly, and to close, we note that the scheme developed in the previous text is suited to both coarse and fine grained parallelism in the sense that each of the decoupled problems can be solved in parallel (the fine graining). Because these problems are independent of each other, they can be solved simultaneously on separate hardware (the coarse graining). Copyright