Possible pole problem in discretization of equations in covariant form

[amrueda]

The DG discretization of equations in covariant form introduced in #47 uses the covariant basis vectors in spherical coordinates to compute the interface fluxes at each RHS evaluation: the state quantities are transformed at the interfaces from each element's local reference frame to the neighbor element's local reference frame passing through the spherical coordinate system.

Since the spherical coordinate system has a singularity at the poles, the transformation of quantities to and from the spherical coordinate system might be subject to instabilities at the poles. Note: We have not observed a case where this transformation is an issue in the current serial implementation. However, I decided to open this issue to investigate the matter in more detail.

Some details:

The transformation at the interfaces is done as follows (citing a comment by @tristanmontoya):

Let $L$ denote the left covariant coordinate system, $R$ the right covariant coordinate system, and $S$ the spherical coordinate system. To transform the vector variables $(v^1, v^2)$ from the left reference frame into the right covariant reference frame we apply the transformation matrix $A_{L \to R}$:

$$(v^1, v^2)_R = A_{L \to R} (v^1, v^2)_L$$

where

$$A_{L \to R} = [A_{R\to S}]^{-1} [A_{L \to S}].$$

The matrices $A_{R\to S}$ and $A_{L \to S}$ contain the covariant basis vectors in spherical coordinates. These matrices are computed from the tree corner node coordinates (v1, v2, v3, v4) and the reference frame coordinates (xi1 and xi2) using the expressions of Guba et al., which are implemented here.

The following example (taken from here) shows that the computation of these matrices is unstable for small perturbations:

If we compute the matrices with

v1 = [0.0, 0.0, 1.0]
v2 = [1.0, 0.0, 0.0]
v3 = [0.0, 1.0, 0.0]
v4 = [0.0, sqrt(2.0) * 0.5, sqrt(2.0) * 0.5]
xi1, xi2 = -1.0, -1.0

then we obtain the covariant basis

julia> basis_covariant
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  0.0   0.353553
 -0.5  -8.96727e-18

However, if we run the code with

v1 = [eps(1.0), eps(1.0), 1.0]
v2 = [1.0, 0.0, 0.0]
v3 = [0.0, 1.0, 0.0]
v4 = [0.0, sqrt(2.0) * 0.5, sqrt(2.0) * 0.5]
xi1, xi2 = -1.0, -1.0

then we get the (very different) matrix

julia> basis_covariant
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 -0.353553   0.25
 -0.353553  -0.25

In other words, small perturbations of the node coordinates on one side of the interface might cause problems in the transformation of the solution.

[OsKnoth]

We can write the transformation as the Jacobian from cartesian to the sphere and than mutiply with a rotation matrix to rotate from Cartesian to spherical.