Skip to content

Differential Geometry#

Problem Statement#

When working with Spectral Element Methods, we are able to handle complex geometry through the use of an unstructured mesh. Additionally, each element in the mesh is allowed to have curved sides. To make this possible, we define a mapping from physical space to a reference "computational" element. The mapping is then approximated by an interpolant whose knots coincide with those used to represent the solution; this is called an isoparameteric mapping.

In SELF, all of the elements are tensor-product element; in 1-D we work with line segments, in 2-D we have quadrilaterals, and in 3-D we have hexahedra. The reference space is then

ξ[1,1]D

where ξ=i=1Dξiξi^ and D is the number of spatial dimensions. The mapping from physical space to computational space is defined as

x=x(ξ)

Recall that we are interested in approximating solutions to the conservation law

st+f=q

The problem we have is that we want to compute the divergence of the conservative flux in physical space, but we only know the coordinate transformation and how to estimate derivatives in computational space. For example, consider calculating sx. Using chain rule, we have

sx=i=1Dsξiξix

The problem here is that we know how to compute xξi, but not the other way around.

These notes explain how we can calculate derivative operators in physical space given that we have the ability to calculate derivatives in computational space and that we know the coordinate transformation from computational to physical space.

Div, Grad, and Curl in mapped coordinates#

In this section, we derive formulas for computing derivatives, divergence, gradient, and curl under coordinate transformations. We then use the formulas to develop algorithms for implementing spectral element methods under coordinate transformations.

We begin with the premise that we have a logically hexahedral space. The physical positions in this space is represented by a vector with three components

x=xjx^jj=1:3

We assume that the physical positions can be calculated in terms of “computational” positions, ξi[1,1]3

x=x(ξ)

Coordinate System Mapping Schematic

When working in mapped coordinates, we want to understand how to calculate divergence, gradient, and curl of scalar and vector quantities. We are usually given a coordinate transformation x=x(ξ) and functions that depend of the computational coordinates ξ.

As an example, to calculate the derivative of a function f, we write

fx=fξ1ξ1x+fξ2ξ2x+fξ3ξ3x=ξfξx

From the above equation, you can show that the divergence of a vector is

f=i=13fξiξi

where

ξi=ξix1x^1+ξix2x^2+ξix3x^3

Similarly, the gradient of a scalar function is

f=i=13fξiξi

and the curl of a vector is

×f=i=13fξi×ξi

The vectors are called the contravariant basis vectors and are usually denoted ai=ξi.

The problem with these formulations is that we usually know the transformation from computational space to physical space ( x=x(ξ) ) and not the other way around. Because of this, it is not readily apparent how to evaluate the contravariant basis vectors. However, we'll show that we can calculate the contravariant basis vectors using quantities we can readily calculate.

Let's first see how we can construct measures of length area and volume. We know that we can relate a small change in x to changes in ξ through the coordinate transformation.

dx=xξ1Δξ1+xξ2Δξ2+xξ3Δξ3=i=13xξiΔξi

The vectors xξi are called the covariant basis vectors and are denoted ai=xξi. This equation can be used to calculate a measure of arc length

ds=|dx|=i=13j=13aiajΔξiΔξj=i=13j=13gi,jΔξiΔξj

The quantities gi,j=aiaj are the elements of the covariant metric tensor.

A small unit area can be calculated by taking the cross product of two units of length.

dAi=(xξjΔξj)×(xξkΔξk)=(aj×ak)ΔξjΔξki,j,kcyclic

A differential volume element is calculated by projecting a differential area in the direction of ai

dV=dAiaiΔξi=(aj×ak)aiΔξiΔξjΔξk=JΔξiΔξjΔξki,j,kcyclic

The quantity (aj×ak)ai is called the Jacobian of the coordinate transformation.

Let’s now look at how to calculate the divergence. The divergence is defined as

f=limΔV01ΔVΔVfdA

The integral on the right hand side is a boundary integral over the faces of the volume element ΔV.

ΔVfdA=(f(a2×a3)|(ξ1+Δξ1,ξ2,ξ3)f(a2×a3)|(ξ1,ξ2,ξ3))Δξ2Δξ3+(f(a3×a1)|(ξ1,ξ2+Δξ2,ξ3)f(a3×a1)|(ξ1,ξ2,ξ3))Δξ3Δξ1+(f(a1×a2)|(ξ1,ξ2,ξ3+Δξ3)f(a1×a2)|(ξ1,ξ2,ξ3))Δξ1Δξ2

With ΔV=JΔξ1Δξ2Δξ3, the divergence becomes

f=1JlimΔV0f(a2×a3)|(ξ1+Δξ1,ξ2,ξ3)f(a2×a3)|(ξ1,ξ2,ξ3)Δξ1+f(a3×a1)|(ξ1,ξ2+Δξ2,ξ3)f(a3×a1)|(ξ1,ξ2,ξ3)Δξ2+f(a1×a2)|(ξ1,ξ2,ξ3+Δξ3)f(a1×a2)|(ξ1,ξ2,ξ3)Δξ3

Upon taking the limit as the volume shrinks to zero, we have

f=1Ji=13[ξi(f(aj×ak))]i,j,kcyclic

Now, we know that the divergence of a constant vector is equal to zero. This implies that

c=0=1Jc[i=13ξi(aj×ak)]i,j,kcyclic

In order for the divergence of a constant vector to be valid for any arbitrary constant, it must be that

i=13ξi(aj×ak)=0i,j,kcyclic

This equation is known as the metric identities. Using the metric identities allows us to write the divergence as

f=1Ji=13fξi(aj×ak)i,j,kcyclic

In comparison with our original formulations for the divergence, we can now relate the contravariant basis vectors to the covariant basis vectors

Jai=aj×aki,j,kcyclic

Similar to the definition of the covariant metric tensor, the contravariant metric tensor is define as

gi,j=aiaj

With this, we can now define the divergence, gradient, and curl under a coordinate transformation.

Operator Conservative Form Non-Conservative Form
Divergence 1Ji=13ξi(Jaif) i=13fξiai
Gradient 1Ji=13ξi(Jaif) i=13fξiai
Curl 1Ji=13ξi(Jai×f) i=13fξi×ai

Free Stream Preservation#

Free stream preservation is a property of a numerical method that refers to the ability of the method to satisfy the metric identities.

With collocation style spectral element methods, we approximate the Jacobian-weighted contravariant basis vectors as a polynomial of degree N. When the coordinate mapping is a polynomial of degree N, the covariant basis vectors are also a polynomial of degree N and the Jacobian weighted contravariant basis vectors are of degree 2N. This introduces an aliasing error that causes a lack of satisfaction of the metric identities.

Kopriva (2006) demonstrated the generation of spurious sound waves in an Euler equations solver when the metric identities are not satisfied discretely. Further, he showed that we can use alternate formulations for the contravariant basis vectors that eliminate the aliasing error through basic vector calculus identities.

First, recall that we can write

u×v=×(vu)=×(uv)

These relations follow from the product rule and the fact that the curl of a gradient is identically zero. From this identity, we can write

ξxm×ξxl=×(xlxm)

Note that ξxm×ξxl=(xmξ2xlξ3xmξ3xlξ2)x^1 (xmξ1xlξ3xmξ3xlξ1)x^2 +(xmξ1xlξ2xmξ2xlξ1)x^3

Further, the contravariant basis vectors can be written

Jai=xξj×xξk=(x2ξjx3ξkx3ξjx2ξk)x^1(x1ξjx3ξkx3ξjx1ξk)x^2+(x1ξjx2ξkx2ξjx1ξk)x^3

If we let i=1,2,3, and i,j,k be cyclic, and let n=1,2,3 with n,m,l cyclic we have that

Jani=x^i(ξxm×ξxl)=x^i(ξ×(xlξxm))=x^i(ξ×(xmξxl))

For example, when n=1, we have that

ξx2×ξx3=(x2ξ2x3ξ3x2ξ3x3ξ2)x^1(x2ξ1x3ξ3x2ξ3x3ξ1)x^2+(x2ξ1x3ξ2x2ξ2x3ξ1)x^3

and when i=1, we have

Ja1=xξ2×xξ3=(x2ξ2x3ξ3x3ξ2x3ξ3)x^1(x1ξ2x3ξ3x3ξ2x1ξ3)x^2+(x1ξ2x2ξ3x2ξ2x1ξ3)x^3

By comparison, we can see that Ja11=x^1(ξx2×ξx3).

Through basic vector calculus, we've shown that we can write the Jacobian weighted contravariant basis vectors as either

Jani=x^i(ξ×(xlξxm))

or

Jani=x^i(ξ×(xmξxl))

Both formulations are guaranteed to be Divergence free and are known as the "conservative forms". The "curl invariant" form is obtained by averaging the two, to get

Jani=12x^iξ×[xmξxlxlξxm]

When calculating the metric terms in a numerical method, we opt to use the curl form of the metric terms, since we can guarantee that the metric terms will be satisfied.