Covariant form as variable-coefficient problem using auxiliary variables

[tristanmontoya]

This potentially replaces #31 by implementing the strategy described in #46, where the geometric information needed by covariant-form fluxes and source terms is passed in as auxiliary variables. The cache is not passed around anymore.

[codecov]

Codecov Report

Attention: Patch coverage is 97.41935% with 8 lines in your changes missing coverage. Please review.

Project coverage is 87.49%. Comparing base (49ec206) to head (7a6146f).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/equations/equations.jl	90.47%	2 Missing ⚠️
src/equations/shallow_water_3d.jl	90.47%	2 Missing ⚠️
...em_p4est/containers_2d_manifold_in_3d_covariant.jl	97.70%	2 Missing ⚠️
...vers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl	97.53%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #47      +/-   ##
==========================================
+ Coverage   84.71%   87.49%   +2.78%     
==========================================
  Files          11       17       +6     
  Lines        1125     1423     +298     
==========================================
+ Hits          953     1245     +292     
- Misses        172      178       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

---

🚨 Try these New Features:

• Flaky Tests Detection - Detect and resolve failed and flaky tests

[amrueda]

Great work, Tristan!
I think that the use of the auxiliary variables makes the implementation cleaner!
Here a couple of comments and questions.

[amrueda]

Does the use of e_lon and e_lat imply a problem in (the vicinity of) the pole singularities? if it does, it should be properly stated. Moreover, is this the standard way to do this? Wouldn't it be a bit more general to use Cartesian unit vectors to store the covariant basis vectors?

[tristanmontoya]

This is the standard way of doing things, and does not introduce any pole problem. This matrix is, in my opinion at least, more convenient than one using Cartesian coordinates, because it converts to/from spherical velocity components, which are typically more meaningful than Cartesian ones in a geophysical context (e.g. initial conditions are commonly defined in terms of zonal and meridional components, not Cartesian ones). It also is 2x2 instead of 3x2 and thus (marginally) cheaper for 2D manifolds.

[amrueda]

You are right, the initial conditions are commonly defined in terms of zonal and meridional components. However my concern is the following:

When I run the code above 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 I get 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 I 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 I get the 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

Very different(!!). Isn't that problematic?... This transformation matrix is used very often to transform the state quantities at the interfaces. What happens if some interface nodes are directly at / in the vicinity of the pole singularity?

[tristanmontoya]

Interesting, I will have to investigate. I see this is used in several other codes, but I don't know if it's something that could cause problems or not.

[tristanmontoya]

If we let $L$, $R$, and $G$ denote the left, right, and global coordinate systems, what ends up happening with the transformation matrices is that you have

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

where

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

If the matrices $A_{R\to G}$ and $A_{L \to G}$ are not (nearly) singular, does it actually matter that there is a pole? The matrices in your example are not close to singular, but perhaps the sensitivity to the node coordinates could still cause issues. The question is what other codes which use a similar element-local mapping (e.g. BRIDGE or CAM/E3SM) do, i.e., is there anything that has to be done to avoid pole issues. The Fortran implementation in E3SM is available here (they call this matrix D instead of A):

https://github.com/E3SM-Project/E3SM/blob/master/components/homme/src/share/cube_mod.F90#L696-L761

[amrueda]

The risk, in my opinion, is that the global coordinate system might be differently defined for two elements at different sides of the interface due to machine precision errors.

[tristanmontoya]

I see. Could it be that even if $A_{L \to R}$ and $A_{R \to L}$ are inverses of each other up to roundoff (this is what we need for conservation) there still may be excessive error introduced that doesn't show up as conservation error? I presume this would be a problem in other codes too.

[tristanmontoya]

I wonder if we should make a separate PR for improvements to this pole issue, as in my opinion, this is more of a scientific problem to be investigated, rather than an implementation error.

[amrueda]

I opened the issue #49
I suggest that we merge this PR but leave this conversation "unresolved" for the moment.

[amrueda]

LGTM! Thank you Tristan!