These html pages are based on the PhD thesis "Cluster-Based Parallelization of Simulations on Dynamically Adaptive Grids and Dynamic Resource Management" by Martin Schreiber.
There is also more information and a PDF version available.

[next] [tail] [up]

A.1 Hyperbolic PDEs

A.1.1 Gauss Lobatto Points

The following table gives examples of basis functions based on Gauss-Lobatto points and their nodal points up to order 2.

Polynomial Nodal point
Degree 0: φ0(x,y) := 1 (1 3,1 3)
Degree 1: φ0(x,y) := 1 - x - y (1 2,1 2)
φ1(x,y) := x (1 2,0)
φ2(x,y) := y (0,1 2)
Degree 2: φ0(x,y) := 1 - 3x + 2x2 - 3y + 4xy + 2y2 (0,0)
φ1(x,y) := 4x - 4x2 - 4xy (1 2,0)
φ2(x,y) := -x + 2x2 (1,0)
φ3(x,y) := 4y - 4xy - 4y2 (0,1 2)
φ4(x,y) := 4xy (1 2,1 2)
φ5(x,y) := -y + 2y2 (0,1)

A.1.2 Jacobi polynomials

The orthogonal Jacobi polynomials on triangle basis up to degree 2 are given in the following Table. Compared to the Gauss-Lobatto Points, the Jacobi polynomials are constructed hierarchically.

Polynomial
Degree 0: φ0(x,y) := √- 2
Degree 1: φ0(x,y) := √- 2
φ1(x,y) := -2 + 6y
φ2(x,y) := 2√- 3(2x - 1 + y)
Degree 2: φ0(x,y) := √- 2
φ1(x,y) := -2 + 6y
φ2(x,y) := (1 - 8y + 10y2)√- 6
φ3(x,y) := 2√- 3(2x - 1 + y)
φ4(x,y) := 3√- 2(-1 + 5y)(2x - 1 + y)
φ5(x,y) := (1 - 2y + y2 - 6x + 6xy + 6x2)√-- 30

A.1.3 Mass Matrix

For nodal basis functions of degree 1 and 2 with Gauss-Lobatto points, the inverse mass matrices are given by

⌊ 72 -3 12 - 3 12 12 ⌋ ||| -3 39 -3 - 27 - 27 12 ||| ⌊ 18 -6 -6 ⌋ ||| 2 4 4 ||| |||⌈ -6 18 -6 |||⌉ ||| 12 -327 72 1239 -327 12 ||| . -6 -6 18 |||| -3 - 4- 12 2- -4- - 3|||| |⌈ 12 - 274- -3 - 274 392 - 3|⌉ 12 12 12 - 3 -3 72

Using orthonormal triangle basis functions based on normalized Jacobi Polynomials, the inverse mass matrix is identical to the identity matrix for arbitrary degree.

A.1.4 Stiffness matrices

For Gauss-Lobatto nodal points, this leads to stiffness matrices each with a single zero-row for degree 1:

⌊ -1∕6 -1∕6 -1∕6⌋ ⌊ -1∕6 -1∕6 -1∕6⌋ ||| 1∕6 1∕6 1∕6 ||| ||| 0 0 0 ||| ⌈ ⌉ ⌈ ⌉ 0 0 0 1∕6 1∕6 1∕6

Computing the matrices x and y for degree 2, this still leads to almost dense matrices.

⌊ ⌋ ⌊ ⌋ || -1∕15 -1∕10 1∕30 - 1∕10 1∕30 1∕30 || || -1∕15 - 1∕10 1∕30 -1∕10 1∕30 1∕30 || ||| 1∕10 0 -1∕10 2∕15 -2∕15 0 ||| ||| 1∕30 -145 -1∕15 -2∕15 - 415- 1∕30 ||| ||| -1∕30 1∕10 1∕15 - 1∕30 1∕10 -1∕30 ||| ||| 0 0 0 0 0 0 ||| . |||| 1∕30 -2∕15 1∕30 -145 - 415 -1∕15 |||| |||| 1∕10 2∕15 0 0 -2∕15 -1∕10 |||| ||⌈ -1∕30 2∕15 -1∕30 415- 415- 1∕15 ||⌉ ||⌈ -1∕30 415- 1∕15 2∕15 145 -1∕30 ||⌉ 0 0 0 0 0 0 -1∕30 - 1∕30 -1∕30 1∕10 1∕10 1∕15

The matrices created for the orthonormal basis function have a higher sparsity pattern. Stiffness matrices for degree 1 are given by

⌊| 0 0 0 ⌋| ⌊| 0 0 0 ⌋| ||⌈ 0 0 0 ||⌉ ||⌈ 3√2- 0 0 ||⌉ 2√3√2 0 0 √3√2 0 0

and for degree 2

⌊ 0 0 0 0 0 0 ⌋ ⌊ 0 0 0 0 0 0 ⌋ |||| 0 0 0 0 0 0 |||| |||| 3√2 0 0 0 0 0 |||| ||| 0 0 0 0 0 0 ||| ||| -4∕3√3 10∕3√6 0 0 0 0 ||| ||| 2√3√2 0 0 0 0 0 ||| ||| √3√2 0 0 0 0 0 ||| . |||| 4 5√2 0 0 0 0 |||| |||| 2 5∕2√2 0 5∕2√3√2 0 0 |||| ⌈ √-- ⌉ ⌈ √-- √-- √-- ⌉ 0 0 0 310 0 0 2∕3 15 -1∕6 30 0 3∕2 10 0 0

Those matrices have clearly a sparser layout and are thus better suited for computation considering the stiffness matrices only so far.
However, using modal basis functions, they have to be transferred to nodal basis functions. Using Gauss-Lobatto nodal points and 1st degree basis functions, this results in the following three matrices:

√ - √- ⌊| 2 -2 -2 3 ⌋| ⌊| 0 0 0 ⌋| ⌊| 0 0 0 ⌋| ||⌈√2- -2 2√3 ||⌉ ||⌈ 0 0 0 ||⌉ ||⌈ 1 1 1 ||⌉ √2- 4 0 2∕3√3 2∕3√3 2∕3√3 1∕3√3 1∕3√3- 1∕3√3

and to the following matrices for 2nd degree

⌊ √2 -2 √6 -2√3 3√2 √30- ⌋ ⌊ 0 0 0 0 0 0 ⌋ ||| √2 -2 √6 0 0 -1∕2√30 ||| ||| 0 0 0 0 0 0 ||| |||| √2 -2 √6 2√3 -3√2- √30- |||| |||| 0 0 0 0 0 0 |||| ||| √2 1 - 1∕2√6 -√3 -9∕4√2 1∕4√30- ||| ||| 0 2∕3√3 0 2∕3√3 2∕3√3 0 ||| ||| √- √- √- √- √-- ||| ||| -1∕4√2 0 -1∕4√2 √2 √2 1∕2√2 ||| |⌈ √2 1 - 1∕√2-6 3 9∕4 2 1∕4 30 |⌉ |⌈ 3√-- 3-√-- √-- √-- |⌉ ⌊ 0 2 4 0 3 6 0 0 0 0 0 0 0 ⌋ - 20 30 0 20 30 -1∕5 30 1∕530 0 . || 0 1 0 1 1 0 || |||| √- √- √- √- |||| ||| -1∕6 6 -2∕3√-6 -1∕6 6 0√- 0√- 1∕3 6 ||| ||| 0√- 1∕3 3 0√- 1∕3 3 1∕3√-3 0√- ||| ||⌈ -1∕2 2 0 1∕42 0 2 1∕4 2 ||⌉ -1∕15√30 2∕15√30 1∕12√30 0 1∕5√30 -610√30

A.1.5 Flux matrices

Only a subset of conserved quantities is involved in the flux computation, see the following matrices for GL nodal points of order 2:

⌊ ⌋ ⌊ ⌋ ⌊ ⌋ | 0 0 1 0 0 0| | 1 0 0 0 0 0 | | 1 0 0 0 0 0 | ||⌈ 0 0 0 0 1 0||⌉ ||⌈ 0 0 0 1 0 0 ||⌉ ||⌈ 0 1 0 0 0 0 ||⌉ . 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0
Thus each row selects the conserved quantities for each nodal point. For modal representation, a similar approach to the stiffness matrices has to be used.

After flux computations, the flux updates only involve the nodal points of the edge for which the flux was computed for. Examples for such matrices generated for flux computations and Gauss Lobatto nodal points in the order of hypotenuse, right and left edge are given by

⌊ 0 0 0 ⌋ ⌊ -2∕15 -1∕15 1∕30 ⌋ ⌊ -2∕15 -1∕15 1∕30 ⌋ ||| 0 0 0 ||| ||| 0 0 0 ||| ||| -1∕15 - 815 -1∕15 ||| ||||-2∕15√2 -1∕15√2 1∕30√2 |||| ||| 0 0 0 ||| ||| 1∕30 -1∕15 -2∕15 ||| ||| 0 0 0 ||| ||| -1∕15 - 8- -1∕15 ||| ||| 0 0 0 ||| . ||| √- 8√ - √-||| |||| 15 |||| |||| |||| |⌈-1∕15√2 -15√2 -1∕15√2|⌉ ⌈ 0 0 0 ⌉ ⌈ 0 0 0 ⌉ 1∕30 2 -1∕15 2 -2∕15 2 1∕30 -1∕15 -2∕15 0 0 0

A.1.6 Butcher tableau

For RK with 2 stages, one possible tableau is

a1,1 := 0 a1,2 := 0
a2,1 := 12 a2,2 := 0
b 1 := 0 b2 := 1

yielding

V0 := U (t) V1 := V0 + Δt (a1,1D1 + a1,2D2 ) = V0 D1 := R (V1) V := V + Δt (a D + a D ) = V + Δt 1D . (A.1) 2 0 2,1 1 2,2 2 0 2 1 D2 := R (V2) U (t+ Δt ) := V0 + Δt (b1D1 + b2D2 ) = V0 + ΔtD2

Another formulation for 2nd order accuracy is given by Heun’s method, yielding the tableau

a1,1 := 0 a1,2 := 0
a2,1 := 1 a2,2 := 0
b1 := 1 2 b2 := 1 2

The tableau for RK3 [But64] yielding 3-rd order accuracy is given by

a1,1 := 0 a1,2 := 0 a1,3 := 0
a2,1 := 1 2 a2,2 := 0 a2,3 := 0
a3,1 := -1 a3,2 := 2 a3,3 := 0
b1 := 1 6 b2 := 2 3 b3 := 1 6

The tableau for classical RK4 yielding 4-th order accuracy [SM03] is given by

a1,1 := 0 a1,2 := 0 a1,3 := 0 a1,4 := 0
a2,1 := 12 a2,2 := 0 a2,3 := 0 a2,4 := 0
a3,1 := 0 a3,2 := 1 2 a3,3 := 0 a3,4 := 0
a4,1 := 0 a4,2 := 0 a4,3 := 1 a4,4 := 0
b1 := 1 6 b2 := 1 3 b3 := 1 3 b4 := 1 6

A.1.7 Rotational invariance of Euler equations

Here, we present the rotational invariancy of the Euler equations [Tor01]: For an edge normal ⃗ne pointing towards (cos(α),sin(α))T , the rotation matrix R(α)e is given by

⌊ 1 0 0 0⌋ | | || 0 cos(α) - sin(α) 0|| || 0 sin (α ) cos(α) 0|| ⌈ ⌉ 0 0 0 1

We test for rotational invariancy by putting the flux terms from Eq. (1.9) into the rotational invariancy formula given in Eq. 2.9 which yields for the left hand side

( cos(α)ru ) ( sin(α)rv ) | ( )| | | || cos(α) ru2 + p || || rusin(α)v || F(U )⋅⃗ne = || r cos(α) uv || + || sin (α) (rv2 + p)|| ( ) ( ) cos(α)(E + p)u sin (α )(E + p)v
(A.2)

and for the right hand side

R(α)-1F(R (α )(U )) = R (α)- 1F((ρ,cos(α)ρu+ sin(α)ρv,- sin(α)ρu+ cos(α )ρv,E )T ) (A.3) ( ) | r(cos(α)u + sin (α) v) | ||ru2cos(α)+ ru sin(α)v + cos(α )p|| = || 2 || (A.4) (r cos (α )uv+ rv sin (α )+ sin(α)p ) (E + p)(cos(α)u+ sin(α)v)
which is equal to F(U) ⋅⃗ne using basic trigonometric calculus.

[next] [front] [up]