[2D manifold in 3D] Option to use global Cartesian coordinates to define exact covariant basis vectors

[tristanmontoya]

As mentioned in #49, there is a possible issue regarding the covariant basis vectors not being well defined at the poles when they are expressed in spherical coordinates. This is not the case when they are defined in Cartesian coordinates. The type parameter GlobalCoordinateSystem has now been added to AbstractCovariantEquations for dispatch purposes, in order to select which coordinate system should be used to construct the exact covariant basis vectors when initializing the auxiliary variables. This can be specified by the user in the constructor for the equations by passing in the keyword argument global_coordinate_system, which can be set to GlobalCartesianCoordinates() (the default) or GlobalSphericalCoordinates().

When a global Cartesian coordinate system is used, the initial_condition_gaussian struct now defines the solid body rotation as $\vec{v}(\vec{x}) = \vec{\omega} \times \vec{x}$, where $\vec{\omega}$ is the angular velocity vector, and computes the analytical solution using Rodrigues' rotation formula, bypassing any use of spherical coordinates whatsoever.

Note: I personally do not see any practical downside to always using global Cartesian coordinates to define the covariant basis. It would simplify the code to just remove the option to use global spherical coordinates. However, the spherical approach appears to be more common in existing codes, and it may be valuable to have it as an option (and to facilitate comparisons of the approaches). If you have an opinion on whether to keep the spherical option or not, please comment on this PR. My current inclination is to keep it for now, and perhaps remove it later on if it becomes cumbersome to continue to support it.

[codecov]

Codecov Report

Attention: Patch coverage is 98.14815% with 2 lines in your changes missing coverage. Please review.

Project coverage is 89.46%. Comparing base (8d2dc5a) to head (e5bb517).
Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/equations/equations.jl	89.47%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #54      +/-   ##
==========================================
+ Coverage   89.41%   89.46%   +0.04%     
==========================================
  Files          20       20              
  Lines        1795     1803       +8     
==========================================
+ Hits         1605     1613       +8     
  Misses        190      190              

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

[amrueda]

Great work!
I just tested convergence for modifications of the spherical advection elixir using the cubed sphere and icosahedral grids. The convergence was broken on the main branch for the icosahedral grid, probably due to #49 (this grid has 5 elements with corners at the north pole and 5 elements with corners at the south pole). The use of global Cartesian coordinates seems to fix this!

Edit(!): It turns out that it was the redefinition of initial_condition_gaussian what fixed the problem. The new form with atan(atan(x[2], x[1])) yields convergence on the main branch as well!

Here are the results I obtain with the main branch:

julia> convergence_test("../examples/elixir_spherical_advection_covariant.jl", 4, cells_per_dimension = (3,3))
####################################################################################################
l2
scalar              vcon1               vcon2               
error     EOC       error     EOC       error     EOC       
7.68e+00  -         0.00e+00  -         0.00e+00  -         
4.66e-01  4.04      0.00e+00  NaN       0.00e+00  NaN       
2.58e-02  4.18      0.00e+00  NaN       0.00e+00  NaN       
1.48e-03  4.12      0.00e+00  NaN       0.00e+00  NaN       

mean      4.11      mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
linf
scalar              vcon1               vcon2               
error     EOC       error     EOC       error     EOC       
6.14e+01  -         0.00e+00  -         0.00e+00  -         
6.55e+00  3.23      0.00e+00  NaN       0.00e+00  NaN       
7.56e-01  3.12      0.00e+00  NaN       0.00e+00  NaN       
4.93e-02  3.94      0.00e+00  NaN       0.00e+00  NaN       

mean      3.43      mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------

julia> convergence_test("../examples/elixir_spherical_advection_covariant_icosahedron.jl", 4, cells_per_dimension = (1,1))
####################################################################################################
l2
scalar              vcon1               vcon2               
error     EOC       error     EOC       error     EOC       
8.74e+00  -         0.00e+00  -         0.00e+00  -         
2.46e+00  1.83      0.00e+00  NaN       0.00e+00  NaN       
1.21e+00  1.02      0.00e+00  NaN       0.00e+00  NaN       
6.53e-01  0.89      0.00e+00  NaN       0.00e+00  NaN       

mean      1.25      mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
linf
scalar              vcon1               vcon2               
error     EOC       error     EOC       error     EOC       
8.08e+01  -         0.00e+00  -         0.00e+00  -         
3.84e+01  1.07      0.00e+00  NaN       0.00e+00  NaN       
2.43e+01  0.66      0.00e+00  NaN       0.00e+00  NaN       
1.62e+01  0.58      0.00e+00  NaN       0.00e+00  NaN       

mean      0.77      mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("scalar", "vcon1", "vcon2")
  :l2        => [1.24716, NaN, NaN]
  :linf      => [0.771628, NaN, NaN]

And here are the results I obtain with tm/new_covariant_metrics:

julia> convergence_test("../examples/elixir_spherical_advection_covariant.jl", 4, cells_per_dimension = (3,3))
####################################################################################################
l2
scalar              vcon1               vcon2               vcon3               
error     EOC       error     EOC       error     EOC       error     EOC       
7.68e+00  -         0.00e+00  -         0.00e+00  -         0.00e+00  -         
4.66e-01  4.04      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
2.58e-02  4.18      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
1.48e-03  4.12      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       

mean      4.11      mean      NaN       mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
linf
scalar              vcon1               vcon2               vcon3               
error     EOC       error     EOC       error     EOC       error     EOC       
6.14e+01  -         0.00e+00  -         0.00e+00  -         0.00e+00  -         
6.55e+00  3.23      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
7.56e-01  3.12      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
4.93e-02  3.94      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       

mean      3.43      mean      NaN       mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("scalar", "vcon1", "vcon2", "vcon3")
  :l2        => [4.11262, NaN, NaN, NaN]
  :linf      => [3.42795, NaN, NaN, NaN]

julia> convergence_test("../examples/elixir_spherical_advection_covariant_icosahedron_newbranch.jl", 4, cells_per_dimension = (1,1))
####################################################################################################
l2
scalar              vcon1               vcon2               vcon3               
error     EOC       error     EOC       error     EOC       error     EOC       
7.49e+00  -         0.00e+00  -         0.00e+00  -         0.00e+00  -         
5.18e-01  3.85      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
4.76e-02  3.45      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
3.67e-03  3.69      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       

mean      3.66      mean      NaN       mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
linf
scalar              vcon1               vcon2               vcon3               
error     EOC       error     EOC       error     EOC       error     EOC       
6.64e+01  -         0.00e+00  -         0.00e+00  -         0.00e+00  -         
1.35e+01  2.29      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
1.67e+00  3.02      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       
1.33e-01  3.66      0.00e+00  NaN       0.00e+00  NaN       0.00e+00  NaN       

mean      2.99      mean      NaN       mean      NaN       mean      NaN       
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("scalar", "vcon1", "vcon2", "vcon3")
  :l2        => [3.66488, NaN, NaN, NaN]
  :linf      => [2.98918, NaN, NaN, NaN]

[amrueda]

Thanks for adding this, @tristanmontoya!

I think that it's good to have the option to choose between GlobalCartesianCoordinates and GlobalSphericalCoordinates.

Here is a preliminary review of the changes.

[amrueda]

LGTM!