Description

Equation
Martin Pham
Fall 2017
Acknowledgements: Keegan Kirk and Abdullah Sivas for time and guidance debugging
1 Formulation
1.1 Computational Domain
For a space-time variable x = (t,x¯) and some constant a, consider the advection problem
∂tu + a∂xu = f E ⊂ R1+1
(1)
u(¯x,t0) = u0(¯x) Ω(t) ∈ R
The space-time domain is E ⊂ R2 and the spatial flow domain Ω at some time t is
Ω(t) = {x¯ ∈ R : (t,x¯) ∈ E}
For a space-time slab En, the bottom and top of the slab are specified by Ω(tn) and Ω(tn+1) respectively. The sides of the time slab are on the boundary of E and are specified by
Qn = ∂E (Ω(tn) ∪ Ω(tn+1))
Let T n = {Kjn} be the tessellation of a flow domain Ω( . We then tesselate a slab En by linear interpolation in time between the spatial tessellations T n and T n+1.
A space-time element Kjn is obtained by connecting the spatial elements Kjn ⊂ Ω(tn) and
1.2 Reference Elements
Let the space-time reference variable be ξ = (ξ0,ξ¯) with space-time reference element Kˆ = (−1,1)2 and spatial reference element Kˆ = (−1,1). At time tn, for some spatial elements Kjn and Kjn+1 we define respectively the mappings

where the shape functions χi are given by

and the points xi(Kjn) and xi(Kjn+1) are given by the movement of the mesh at time tn and tn+1 respectively.
For a space-time element Kjn we define the mapping

1.3 Weak Form
Define the space-time gradient operator and flux
∇ = (∂t,∂x) F(u) = (u,au)
Then the differential equation in (1) becomes
∇ · F = f (2)
(3)
The second integral term is replaced by an approximation of the flux F using a numerical flux Fˆ that will consider both sides of the solution across discontinuities between elements. On a space-time element, define the grid-dependent space-time normal that is consistent with a causality in time argument as

(−1,0)

n = (1,0)
(−vgn,¯ n¯) Kjn(t+n )
Kjn+1(t−n+1)
Qnj
(4)
Since the solution is double-valued at points of discontinuity across elements, we introduce the approximation of u at these points as

Z Z Z Z
vF · nds = − vudˆ x¯ + vudˆ x¯ + v(au− vgu)¯nds
∂K Kjn Kjn+1 Qnj
(5)
Z Z Z
= − vuRdx¯ + vuLdx¯ + v(au− vgu)¯nds
Kjn Kjn+1 Qnj
Denote the jump and average operators as

=v gˆ(uL,uR)ds (6)
J K
Z
= X (vL − vR)¯nLgˆ(uL,uR)ds
S∈Sn S Z
= X (vL − vR)fˆ(uL,uR)ds
S∈Sn S
where fˆ(uL,uR) = gˆ(uL,uR)¯nL. We choose the upwind flux
K
Note that if there is no grid motion, ¯n becomes a unit normal

Note that the left-most facet and right-most facet require special treatment dependant on the direction of the advection, no jump in test functions is imposed.
Consider the problem (3), applying our integral terms from (5) and (6) we obtain

Splitting the computational domain into space-time slabs

(7)
1.4 Polynomial Expansion Coefficients
u,v ∈ V k = {u ∈ L2(Kj) : u|Kj ∈ Pk(Kj)}
for k the highest degree Legendre polynomial basis function on element Kj. We choose a linear polynomial approximation k = 1. On the reference element Kˆj = (−1,1)2, we define the basis functions Ψm as
Ψˆ1 = 1 Ψˆ2 = ξ¯= ξ1 Ψˆ3 = ξ0
where Ψˆ1 determines the element average, Ψˆ2 the slopes in the spatial dimension and Ψˆ3 the slopes in the time dimension. The basis functions on a space-time element Kj are thus

Express the approximations u,v as
3 u(t,x¯)|K = X uˆmΨm(t,x¯) = uˆ1Ψ1 + uˆ2Ψ2 + uˆ3Ψ3
m=1
(8)
3 v(t,x¯)|K = X vˆmΨm(t,x¯) = vˆ1Ψ1 + vˆ2Ψ2 + vˆ3Ψ3 m=1
In order to compute the integrals in the weak form of the problem we must obtain an expression for the derivatives with respect to space-time variables t,x¯ in terms of reference variables ξ0,ξ1. Consider, by chain rule
(9)
Recall that ∆tn = tn+1 − tn and . From the definition of the mapping

The matrix vector equation in (9) becomes

where the matrix J is the Jacobian of the mapping GK for an element Kj. Inverting the equation, we obtain

For every element we apply the relationship
Z Z
f(x)dtdx¯ = f(GnK(ξ))|J|dξ0dξ¯
K Kˆ
(11)

(12)
and the second integral
(13)
For a d dimensional space-time problem there are three d − 1 dimensional integrals in d dimensions in (7)
ZZ
vLuLdx¯(vL − vR)fˆ(uL,uR)ds
Kjn+1S
The first d − 1 integral is an integral over the top of the space-time slab at time tn+1, the spatial element Kjn+1. Let |Kjn+1| be the length of the spatial element. Note that
t = tn+1 =⇒ ξ0 = 1
and so

Taking the derivative with respect to ξ¯, we obtain the Jacobian of the spatial reference mapping

(14)
where because we are at the top of our space-time element at time tn+1 (and so (ξ0)L = 1) we have made the choice
uL = uˆn1 + uˆn2ξ¯+ uˆn3(ξ0)L = uˆn1 + uˆ2nξ¯+ uˆn3
vL = vˆ1n + vˆ2nξ¯+ vˆ3n(ξ0)L = vˆ1n + vˆ2nξ¯+ vˆ3n
The second d−1 integral is an integral over the bottom of the slab at time tn. We apply a similar method as above however because we are at the bottom of our space-time element (and so (ξ0)L = −1), we upwind information from the previous space-time slab evaluating at the corresponding reference time (i.e. (ξ0)R = 1)
uR = uˆ1n−1 + uˆn2−1ξ¯+ uˆn3−1(ξ0)R = uˆn1−1 + uˆn2−1ξ¯+ uˆn3−1
vL = vˆ1n + vˆ2nξ¯+ vˆ3n(ξ0)L = vˆ1n + vˆ2nξ¯− vˆ3n
(15)
The third d − 1 integral is an integral over the shared facet of elements within a slab. In the case of one spatial dimension, the space-time element Kjn has two facets S. Let Sj be the facet corresponding to the discontinuity along the interpolation of point (xnj )± to
. Let ˆvm,j denote the test function expansion coefficients on the jth element. The polynomial expansions for an integral over Sj is then
vL = vˆ1 + vˆ2(ξ¯)L + vˆ3ξ0 = vˆ1,j − vˆ2,j + vˆ3,jξ0 vR = vˆ1 + vˆ2(ξ¯)R + vˆ3ξ0 = vˆ1,j−1 + vˆ2,j−1 + vˆ3,j−1ξ0
Let Sj+1 be the facet corresponding to the discontinuity along the interpolation of point to ( . The polynomial expansions for an integral over Sj+1 is then
vL = vˆ1 + vˆ2(ξ¯)L + vˆ3ξ0 = vˆ1,j + vˆ2,j + vˆ3,jξ0 vR = vˆ1 + vˆ2(ξ¯)R + vˆ3ξ0 = vˆ1,j+1 − vˆ2,j+1 + vˆ3,j+1ξ0
Z Z Z
(vL − vR)fˆ(u)ds = vLfˆ(uL,uR)ds − vRfˆ(uL,uR)ds
Sj Sj Sj
(16)
Z Z Z
(vL − vR)fˆ(u)ds = vLfˆ(uL,uR)ds − vRfˆ(uL,uR)ds
Sj+1 Sj+1 Sj+1
where

(17)
(18)
2 Numerical Results
Consider the manufactured solution
u(¯x,t) = sin(2πx¯)sin(3πt)
over time t ∈ [t0,TF ] = [0,1] with flow domain ¯ ] for the problem
∂tu + ∂xu = 3π sin(2πx)cos(3πt) + 2π cos(2πx)sin(3πt) u(¯x,t0) = u0(¯x) = 0 We consider the L2 error computed by quadrature

and determine orders of convergence using different numbers of space-time elements
Nel = 2,4,8,…,1024
We set the number of time slabs equal to the number of space-time elements per slab. Consider the rates of convergence below for the solution at final time. As the computational domain is refined, second-order accuracy is obtained. This supports the theoretical rates of convergence found in the literature.
Nel L2 error Order
2 0.063742775144965
4 0.071522780207995 -0.166140986989625
8 0.010296406509315 2.796261899135431
16 0.002043067685942 2.333332012569084
32 0.000467315293649 2.128268843104989
64 0.000114501523735 2.029029455315673
128 0.000028575779637 2.002502031855937
256 0.000007153817545 1.998007631687136
512 0.000001790713916 1.998178452160747
1024 0.000000448025504 1.998882104185325

Figure 1: Note that exact solution is plotted in black, the approximate in red