WIP: p4est solver on a spherical shell

[amrueda]

This PR implements a version of the 2D p4est solver that runs on a spherical shell that is constructed from a cubed sphere. Since there are three Cartesian space dimensions on the surface of the sphere, we use the 2D p4est solver in combination with a 3D equation.

The connectivities are inherited from the p8est 3D cubed-sphere implementation, but they are adapted to p4est (2D). The metric terms are computed with a cross-product form, as it is standard for the DG implementations of the shallow-water equations on the surface of a sphere. See, e.g.,

• Giraldo (2001). A spectral element shallow water model on spherical geodesic grids. http://faculty.nps.edu/fxgirald/Homepage/Publications_files/Giraldo_IJNMF_2001.pdf.

I did not see any particular advantages of using a curl form.

The following points need to be addressed (more information in the threads below):

• The current implementation of get_contravariant_vector causes too many allocations.
• Improve call to new source term that depends on du.

The following visualization uses this version of Trixi2Vtk:

earth.mp4

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[amrueda]

I removed this check because the new implementation requires a 2D mesh (P4estMesh{2}) with a 3D equation. Should I just remove this line or maybe change the check?

[ranocha]

That's a good question. I think we should analyze where we use this kind of assumption and try to get it right. For example, I guess everything will fail badly if you use an UnstructuredMesh2D with a 3D equation, won't it?

[amrueda]

Yes, it will fail. We would need special constructors that create those meshes using nodes with three (x, y, z) node coordinates, and we would then need to adjust the metric term computation of all mesh types. It is doable but I think it is out of the scope of the present PR and probably no one will ever use it 🤣.

Then maybe I just check for the particular case of a P4estMesh2D with 3D node coordinates in combination with a 3D equation, allow that combination and bring back the @assert for any other case. Does that sound like a good solution?

[amrueda]

This new version of get_contravariant_vector causes too many allocations. Unfortunately, ndims(contravariant_vectors) - 3 does not work for the simulations on the spherical shell and size() is evaluated in runtime. Is there a better way to do this?

[sloede]

Is this code only called from your code or also from existing routines? In the current implementation it is not type stable, thus the allocations.

We might get around it easily if it is only used from code you control, and harder (but not impossible) if it is in generic routines as well.

[amrueda]

No, this is called in the existing volume and surface integral routines. In fact, I temporarily changed the existing volume integral calls, such that the code runs a bit faster. See, e.g.,

Trixi.jl/src/solvers/dgsem_structured/dg_2d.jl

Lines 78 to 81 in b89e333

	#Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
	Ja11 = contravariant_vectors[1, 1, i, j, element]
	Ja12 = contravariant_vectors[2, 1, i, j, element]
	contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2

[sloede]

Hm. I wonder if what you are doing in this PR is fundamental enough of a conceptual difference to require a new mesh type. That is, to not try to shoehorn this into P4estMesh{2} but rather something like P4estMesh{3/2}?

[ranocha]

Alternatively, we may need to pass an additional argument to get_contravariant_vector that determines the number of dimensions at compile time - the mesh?
This would be a kind of semi-breaking change since it's not documented anywhere as far as I know but used in the internals. Thus, it may break the workflow for people using Trixi.jl as a library.

[sloede]

True. The question is whether the mesh truly "owns" the dimensionality of a problem. But it would certainly be less intrusive than adding a new mesh type.

[ranocha]

Well, it looks like there isn't anything like "the dimensionality of a problem" anymore. We need to distinguish between the dimensionality of the geometry and the equations. From my point of view, the contravariant vectors are clearly associated with the geometry, i.e., the mesh.

However, this is really a tough problem since we want to pass the contravariant vectors as normal directions to the fluxes. Thus, we assume that the dimensionality of the equations and the mesh are the same. This requires some really careful design to achieve a consistent splitting.

[amrueda]

We need to distinguish between the dimensionality of the geometry and the equations.

Exactly!
The contravariant vectors are actually associated with both the dimensionality of the geometry and the equations. For instance:

• For a 2D geometry and a 2D equation, the contravariant vectors are of size 2 x 2 x nnodes x nnodes x nelements.
• For a 2D geometry and a 3D equation, the contravariant vectors are of size 3 x 3 x nnodes x nnodes x nelements.
• For a 3D geometry and a 3D equation, the contravariant vectors are of size 3 x 3 x nnodes x nnodes x nnodes x nelements.

For the first and last cases, ndims(contravariant_vectors) - 3 works. For the (new) second case, it does not.

For this particular function, we would actually need to pass the number of dimensions of the equation to define the size of the SVector... Or maybe there is a workaround. 🤔

[amrueda]

Hm. I wonder if what you are doing in this PR is fundamental enough of a conceptual difference to require a new mesh type. That is, to not try to shoehorn this into P4estMesh{2} but rather something like P4estMesh{3/2}?

This function is mesh type agnostic, so that won't fix this particular problem.

[amrueda]

In amrueda#12, I propose a workaround to this problem using a new type that extends Array. What do you think of this solution? @sloede @ranocha

[amrueda]

Convergence test with examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl. It computes the advection of a Gaussian around the sphere.

• polydeg = 3:

####################################################################################################
l2
rho                 rho_v1              rho_v2              rho_v3              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       
8.45e-03  -         8.23e-03  -         1.84e-03  -         0.00e+00  -         6.44e-02  -         
8.47e-04  3.32      8.17e-04  3.33      2.04e-04  3.17      0.00e+00  NaN       6.46e-02  -0.00     
2.30e-05  5.20      2.26e-05  5.17      4.08e-06  5.65      0.00e+00  NaN       6.46e-02  0.00      
1.06e-06  4.44      9.43e-07  4.59      3.73e-07  3.45      0.00e+00  NaN       6.46e-02  0.00      

mean      4.32      mean      4.36      mean      4.09      mean      NaN       mean      -0.00     
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_v2              rho_v3              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       
9.49e-02  -         9.65e-02  -         1.37e-02  -         0.00e+00  -         4.09e-01  -         
8.23e-03  3.53      8.23e-03  3.55      1.95e-03  2.82      0.00e+00  NaN       4.22e-01  -0.04     
5.89e-04  3.80      5.93e-04  3.79      6.77e-05  4.85      0.00e+00  NaN       4.17e-01  0.02      
3.44e-05  4.10      3.47e-05  4.10      9.36e-06  2.85      0.00e+00  NaN       4.16e-01  0.00      

mean      3.81      mean      3.81      mean      3.51      mean      NaN       mean      -0.01     
----------------------------------------------------------------------------------------------------

• polydeg = 4:

####################################################################################################
l2
rho                 rho_v1              rho_v2              rho_v3              rho_e
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC
2.52e-03  -         2.37e-03  -         8.18e-04  -         0.00e+00  -         6.47e-02  -
3.14e-05  6.33      3.10e-05  6.26      5.44e-06  7.23      0.00e+00  NaN       6.46e-02  0.00
8.98e-07  5.13      8.81e-07  5.14      1.57e-07  5.11      0.00e+00  NaN       6.46e-02  -0.00
5.24e-08  4.10      5.02e-08  4.13      9.30e-09  4.08      0.00e+00  NaN       6.46e-02  0.00

mean      5.19      mean      5.18      mean      5.47      mean      NaN       mean      0.00
----------------------------------------------------------------------------------------------------
linf   
rho                 rho_v1              rho_v2              rho_v3              rho_e
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC
2.04e-02  -         1.99e-02  -         6.87e-03  -         0.00e+00  -         4.20e-01  -
9.62e-04  4.41      9.69e-04  4.36      1.13e-04  5.92      0.00e+00  NaN       4.18e-01  0.01
2.86e-05  5.07      2.87e-05  5.08      4.42e-06  4.68      0.00e+00  NaN       4.16e-01  0.01
1.55e-06  4.20      1.55e-06  4.21      1.66e-07  4.73      0.00e+00  NaN       4.16e-01  0.00

mean      4.56      mean      4.55      mean      5.11      mean      NaN       mean      0.00
----------------------------------------------------------------------------------------------------