Implement a new mesh type that uses p4est to manage cell connectivity and mesh refinement

[efaulhaber]

This issue is a to-do list of my project to implement a new mesh type P4estMesh.

• Implement the new mesh type and corresponding data structures. Create a constructor to build a structured mesh that can be curved by specifying a transformation mapping (just as with CurvedMesh). Only work with periodic boundaries and conforming (uniformly refined) meshes (Implement P4estMesh with periodic boundaries in 2D #578).
• Make SaveSolution and SaveRestart callbacks work with the new mesh type (Implement save/restart with P4estMesh #596).
• Add visualization code for the new mesh type to Trixi2Vtk (Add support for P4estMesh Trixi2Vtk.jl#33).
• Implement non-periodic boundaries. If this doesn't turn out to be more complicated than I think right now, I will also add another mesh constructor that reads Abaqus mesh files. p4est can read these files, and hopefully I can extract the cell coordinates from the generated connectivity easily. This would allow me to test the new mesh type with unstructured meshes, which has already been implemented but never tested. Unfortunately, p4est can't read curved mesh files, so this will be a lot more work to implement, which I won't do in this project (Implement non-periodic boundaries with P4estMesh #598).
• Add support for non-conforming meshes (conforming tree structure but not uniformly refined) (Add support for static non-conforming P4estMeshes (no AMR) #621).
• Implement AMR with p4est (Implement AMR with p4est #618).
• Optimize performance of AMR with P4estMesh (Improve performance of AMR with P4estMesh #627).
• Implement everything except FSP on non-conforming meshes in 3D (Implement 3D support for P4estMesh #637).
• Fix FSP on non-conforming meshes in 3D (Add free stream preservation tests for non-conforming P4estMesh and EC test #666).
• Build the Cubed Sphere mesh with the new mesh type (Implement the Cubed Sphere as P4estMesh #689).
• Add brief docs (add brief docs on P4estMesh #699)

There are a few things to do besides this linear implementation plan:

• Compare performance to the other mesh types.
• Decide whether P4est.jl should be a hard dependency or via Requires.jl (Implement P4estMesh with periodic boundaries in 2D #578 (comment)).
• Decide whether to using or import P4est because Trixi's namespace will be populated with a lot of symbols (Implement P4estMesh with periodic boundaries in 2D #578 (comment)).
• Decide how to proceed with 3D code because most of this will be very similar but with p8est_ instead of p4est_ functions (Implement P4estMesh with periodic boundaries in 2D #578 (comment)).

[sloede]

Thanks a lot for opening this!

[efaulhaber]

As requested by @ranocha in #578, I did a performance analysis.
The benchmark with PkgBenchmark still hadn't finished after I was done mowing the lawn, so I just did a manual benchmark.
These are my observed minimum times of rhs!:

	TreeMesh	CurvedMesh	P4estMesh
elixir_advection_basic	64.300 μs	72.800 μs	101.100 μs
elixir_euler_source_terms	261.200 μs	246.400 μs	378.600 μs

Running elixir_advection_basic until t=100 leads to these timer outputs.
TreeMesh:

 ───────────────────────────────────────────────────────────────────────────────
            Trixi.jl                    Time                   Allocations      
                                ──────────────────────   ───────────────────────
        Tot / % measured:            1.62s / 96.9%           22.2MiB / 75.6%    

 Section                ncalls     time   %tot     avg     alloc   %tot      avg
 ───────────────────────────────────────────────────────────────────────────────
 rhs!                    20.0k    1.34s  85.6%  67.1μs   12.2MiB  72.8%     641B
   volume integral       20.0k    746ms  47.5%  37.3μs     0.00B  0.00%    0.00B
   interface flux        20.0k    226ms  14.4%  11.3μs     0.00B  0.00%    0.00B
   surface integral      20.0k    147ms  9.39%  7.37μs     0.00B  0.00%    0.00B
   prolong2interfaces    20.0k    128ms  8.13%  6.38μs     0.00B  0.00%    0.00B
   Jacobian              20.0k   45.6ms  2.91%  2.28μs     0.00B  0.00%    0.00B
   ~rhs!~                20.0k   30.4ms  1.94%  1.52μs   5.81MiB  34.6%     305B
   reset ∂u/∂t           20.0k   8.41ms  0.54%   421ns     0.00B  0.00%    0.00B
   prolong2mortars       20.0k   4.24ms  0.27%   212ns   6.41MiB  38.2%     336B
   prolong2boundaries    20.0k   2.55ms  0.16%   127ns     0.00B  0.00%    0.00B
   mortar flux           20.0k   2.20ms  0.14%   110ns     0.00B  0.00%    0.00B
   boundary flux         20.0k   1.84ms  0.12%  91.8ns     0.00B  0.00%    0.00B
   source terms          20.0k   1.13ms  0.07%  56.3ns     0.00B  0.00%    0.00B
 I/O                        42    186ms  11.9%  4.44ms   4.13MiB  24.6%   101KiB
 analyze solution           41   38.5ms  2.45%   940μs    444KiB  2.58%  10.8KiB
 calculate dt            4.00k   1.54ms  0.10%   384ns     0.00B  0.00%    0.00B
 ───────────────────────────────────────────────────────────────────────────────

CurvedMesh:

 ─────────────────────────────────────────────────────────────────────────────
           Trixi.jl                   Time                   Allocations      
                              ──────────────────────   ───────────────────────
       Tot / % measured:           1.86s / 97.2%           17.8MiB / 69.6%    

 Section              ncalls     time   %tot     avg     alloc   %tot      avg
 ─────────────────────────────────────────────────────────────────────────────
 rhs!                  20.0k    1.53s  84.5%  76.4μs   5.80MiB  46.8%     304B
   volume integral     20.0k    944ms  52.3%  47.2μs     0.00B  0.00%    0.00B
   interface flux      20.0k    319ms  17.7%  15.9μs     0.00B  0.00%    0.00B
   surface integral    20.0k    154ms  8.53%  7.70μs     0.00B  0.00%    0.00B
   Jacobian            20.0k   74.0ms  4.10%  3.70μs     0.00B  0.00%    0.00B
   ~rhs!~              20.0k   23.8ms  1.32%  1.19μs   5.80MiB  46.8%     304B
   reset ∂u/∂t         20.0k   8.37ms  0.46%   419ns     0.00B  0.00%    0.00B
   boundary flux       20.0k   2.08ms  0.12%   104ns     0.00B  0.00%    0.00B
   source terms        20.0k   1.28ms  0.07%  64.1ns     0.00B  0.00%    0.00B
 I/O                      42    176ms  9.73%  4.19ms   4.11MiB  33.1%   100KiB
 calculate dt          4.00k   55.7ms  3.08%  13.9μs     0.00B  0.00%    0.00B
 analyze solution         41   48.3ms  2.67%  1.18ms   2.50MiB  20.2%  62.5KiB
 ─────────────────────────────────────────────────────────────────────────────

P4estMesh:

 ───────────────────────────────────────────────────────────────────────────────
            Trixi.jl                    Time                   Allocations      
                                ──────────────────────   ───────────────────────
        Tot / % measured:            2.28s / 97.7%           14.0MiB / 61.4%    

 Section                ncalls     time   %tot     avg     alloc   %tot      avg
 ───────────────────────────────────────────────────────────────────────────────
 rhs!                    20.0k    2.14s  95.7%   107μs   6.11MiB  70.8%     320B
   volume integral       20.0k    957ms  42.9%  47.8μs     0.00B  0.00%    0.00B
   interface flux        20.0k    584ms  26.2%  29.2μs     0.00B  0.00%    0.00B
   prolong2interfaces    20.0k    319ms  14.3%  15.9μs     0.00B  0.00%    0.00B
   surface integral      20.0k    156ms  7.00%  7.81μs     0.00B  0.00%    0.00B
   Jacobian              20.0k   74.8ms  3.35%  3.74μs     0.00B  0.00%    0.00B
   ~rhs!~                20.0k   35.7ms  1.60%  1.79μs   6.11MiB  70.8%     320B
   reset ∂u/∂t           20.0k   8.58ms  0.38%   429ns     0.00B  0.00%    0.00B
   source terms          20.0k   1.29ms  0.06%  64.4ns     0.00B  0.00%    0.00B
 calculate dt            4.00k   55.3ms  2.48%  13.8μs     0.00B  0.00%    0.00B
 analyze solution           41   41.1ms  1.84%  1.00ms   2.52MiB  29.2%  62.9KiB
 ───────────────────────────────────────────────────────────────────────────────

The volume integral is the same as in CurvedMesh.
Interface flux and prolong2interfaces is slower than in both TreeMesh and CurvedMesh, which should be investigated.

[ranocha]

Thanks for running these benchmarks! It's a good smoke test for now and nothing completely suspicious popped up in my opinion. We can keep the performance differences in mind and investigate them later.

[sloede]

I think the most interesting comparison will be against the unstructured solver, since clearly the difference in nodal orientation has an non-negligible impact on performance.

But I agree, nothing to tackle immediately but something to keep in mind.

[ranocha]

Anything left to do or can we close this, @efaulhaber?

[efaulhaber]

Yes, the new mesh type has been implemented.
I removed the point "Do some nice computations on the cubed sphere", as this is not really a part of "implementing a new mesh type".