Consider the system (2.17), expressed as
∂u
∂t +∂Au
∂x +Bu= 0 (3.1)
whereu= (ρ0, u0, p0)T being the primitive variables. This is a system of coupled advection equations with constant matricesAandB, given as
A=
u¯ ρ¯ 0 0 u¯ 1/¯ρ 0 γp¯ u¯
, B=
0 0 0
¯
u¯ux/¯ρ 2¯ux −ρ¯x/¯ρ2 0 (1 +γ)¯px (1 +γ)¯ux
(3.2)
where the matrixB ≡0, if the mean flow is uniform on either side of the heat source. In the following presentation, such a uniform flow is assumed, because it corresponds to the main situation in the acoustic model, and because it does not have change the derivation of the numerical scheme. Hence, the equation (3.1) is presented in a general conservation form
∂u
∂t +∂f(u)
∂x = 0 (3.3)
withf(u) =Au
In the following presentation of the scheme, the equation to be discretized is a scalar equation. This is only to make the notation easier, it is directly applicable on each equation in the system (3.3).
So the problem equation is
∂u
∂t +∂f(u)
∂x = 0 (3.4)
along with some initial conditionu(x,0) =u0(x) and boundary conditions.
3.1 Discontinuous Galerkin Formulation 31
The solution is sought in the spatial domain Ω : x ∈ [0, L], initially by a fi-nite element formulation. The fifi-nite element formulation consists of dividing the spatial domain in non-overlapping and not necessarily uniform elements Dk =
xk−, xk+
, such that Ω =SK
k=1Dk, and expressing the solutionuon each element as a local polynomial approximation ukh. TheN’th order polynomial representation can be either on a modal or a nodal basis, such that the local element-wise solution forx∈Dk is represented as
ukh(x, t) = nodal values. The basis considered here, of which different choices could be made, are the Legendre polynomials ψn and the Lagrange interpolating poly-nomials li(x). The order of the approximation is one of the parameters to control the discretization error, the other being the number of elements cover-ing the domain. This polynomial approximation is what gives the scheme the nice property of having spectral convergence, i.e. for polynomial projections and smooth solutionsu∈Hpwithplarge – which denotes solutions withp−1 continuous derivatives and bounded p’th derivative – the approximation error ku−uhkisO(N−p), which is a very fast convergence.
The solution to the equation is then the sum of these local solutions u(x, t)≈uh(x, t) =
XK k=1
ukh(x, t) . (3.6)
The ”approximation error” by this polynomial representation is the quantity, that is wanted minimized, or more precisely, wanted to be a zero-function, in the weak sense of a L2-norm (the inner product). This is done by expressing a residual R(uh) by the polynomial approximation inserted into the equation, then multiplying with a test function φ, locally defined on each element and continuous on this, and then taking the inner product. In the Galerkin frame-work, these test functions are chosen from the same function space as the one used to express the polynomial approximation (the trial functions).
So on each element, the problem is defined as
∀k, l
is a general representation used henceforth, where the test functionφcould be either the Legendre basis functionsψor the Lagrange basisl, as given in (3.5).
Performing integration by parts on the spatial derivative once, gives the weak form of the scheme
∀k, l dG-method separates itself from the finite element formulation. In the finite element formulation, the basis usually vanishes at element edges, such that the right hand side in (3.9) is zero.
Since the trial space, onto which the numerical solution is decomposed, has no requirement of continuity across elements, this flux can be multiply defined as bothf(uh(xk+, t)) andf(uh(xk+1− , t)), sincexk+=xk+1− . This ambiguity is the reason for introducing the so-callednumerical flux or trace f∗ in (3.9), which is a key concept in the dG-method. How this is defined and what function it serves, will be adressed later. But from the above expression it is apparent, that something having to do with an element edge, comes into the scheme by this numerical flux.
This deficiency of continuity requirement in-between elements is also a reason why the dG method has efficiency advantages. The elements decouple in a sense, so that for large problems (in terms of grid), each element can be solved for, on its own i.e. the scheme is parallizable, a valuable property when dealing with large time-dependent problems.
The dG scheme has a second form, called the strong form, which results from the weak form (3.9) by re-doing integration by parts
∀k, l These two forms are mathematically the same, although the strong form should have a tendency of better convergence. Also it is the form used in the imple-mentation, although no difference was found between using the weak and the strong form.
Inserting the polynomial representation (3.8), and of course a similar representa-tion of the flux funcrepresenta-tion, into the weak form (3.9), and interchanging summarepresenta-tion and integration yields for each element
XN
3.1 Discontinuous Galerkin Formulation 33
or equivalently expressed in matrix-vector notation, with ˆukh = ˆuk0, . . . ,uˆkNT
being the vector of solution coefficients, and the same for ˆfkh: Mkdˆukh
dt − SkTfˆkh=−[f∗φl]x
k +
xk− (3.12)
where the matrices are being defined as Mkij =
The strong form results in a similar form, as Mkduˆkh
dt + Skfˆkh=
f(ukh)−f∗ φlxk+
xk− (3.14)
The matrices Mk and Sk are local operators, since they depends on the basis defined on each element and in principle should be defined for each element.
Without going into details with the implementation, this is overcome by defining the basis on a reference element, such that these matrix operators operate on the reference element and then mapped onto the physical, spatial element. This is however more of an efficiency concern.
The independency of each element equation system – (3.9) or (3.10) – with respect to each other, adds the property of being hp-adaptive to the scheme.
This means that adjacent local solutions can be approximated in different order and/or over elements of different lengths. So the grid spacing and approximation order can be adjusted to parts of the problem, if needed. This can also be done during the calculations, with an extra computational effort of course, so that parts of the solution can be discretized ”at a smaller scale” or eventually tracked.
This above form are reduced even further, again without considering the details about construction etc. such that e.g. for the strong form
dˆukh
dt + Dkfˆkh= Mk−1
f(ukh)−f∗ φlxk+
xk− (3.15)
where the matrix Dk acts as a differentiation matrix (or a interpolation form of the derivative), in much the same way as in spectral methods.
This is the sought semi-discretization of the PDE, now we can apply an integrator for the time-derivative, where the right hand side of this ”ODE” is the spatial discretization, i.e.
dukh
The integrator used is the 4’th order, 5 stage Low Storage Runge-Kutta method presented in [2]. Alternatives have been tried in form of a 4’th order, 6 stage Low Dispersion Low Dissipation Runge-Kutta (also in low storage form), but without distinct differences. So the cheaper 5 stage scheme was used.
The discretization in space sets a bound on the time-step, in a CFL-sense.
The time-stepping should ensure that the domain of dependence (regarding the hyperbolicity of the problem) is respected. So a bound on time-step involves the minimum spatial discretization length, but it also has to take into account the speed of the propagating waves in the system. A bound can be formulated as ([2, p.62]),
∆t≤Cmin ∆x
ρ(A) (3.18)
with ρ(A) being the (global) spectral radius of the flux coefficient matrix in (3.1). When frequency analysis is performed on recorded series, the time-step (or sampling rate) is small enough to guarantee that waves in the frequency range considered are represented.