WIP: modify calc_node_coordinates! to allow for P4estMesh{2} with 3D coordinates

[tristanmontoya]

In TrixiAtmo.jl, we solve PDEs on two-dimensional surfaces in three-dimensional space by using P4estMesh{2} with three-dimensional node coordinates and redefining init_elements accordingly to set up the appropriate containers for such meshes. We currently replace calc_node_coordinates! with a function which is otherwise identical to the corresponding Trixi.jl function, but changes the line

tmp1 = StrideArray(undef, real(mesh),
                       StaticInt(2),
                       static_length(nodes), static_length(mesh.nodes))

to the following:

tmp1 = StrideArray(undef, real(mesh),
                       StaticInt(size(mesh.tree_node_coordinates, 1)),
                       static_length(nodes), static_length(mesh.nodes))

If this doesn't cause any issues elsewhere, it would be helpful for our purposes to have this change be part of Trixi.jl, not only because it avoids repetition of the entire function in TrixiAtmo.jl, but also because Trixi2Vtk.jl calls Trixi.calc_node_coordinates! and thus requires such a change in order for this PR to be used for plotting 2D surfaces in 3D. Similar modifications could be made to other mesh types/dimensions, but since I don't see a use case for that currently, I haven't done that in this PR.

[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 with its PR number.

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.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.30%. Comparing base (afd8060) to head (433fd18).
Report is 10 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2046   +/-   ##
=======================================
  Coverage   96.30%   96.30%           
=======================================
  Files         466      466           
  Lines       37227    37227           
=======================================
  Hits        35850    35850           
  Misses       1377     1377           

Flag	Coverage Δ
unittests	96.30% <ø> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

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

[ranocha]

Did you check the performance and type stability of this?

[tristanmontoya]

I hadn't thought about the fact that this would get called repeatedly when using mesh adaptation, so I initially wasn't too concerned about performance. I've now done some tests with a 2D mesh.

julia> using Trixi, BenchmarkTools

julia> using StaticArrayInterface: static_length

julia> using StrideArrays: StrideArray, StaticInt

julia> solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs);

julia> mesh = P4estMesh((8, 8), polydeg = 3,
                        coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0),
                        initial_refinement_level = 1);

Existing code:

julia> @benchmark StrideArray(undef, real(mesh),
                              StaticInt(2),
                              static_length(solver.basis.nodes), static_length(mesh.nodes))
BenchmarkTools.Trial: 10000 samples with 195 evaluations.
 Range (min … max):  489.103 ns … 539.166 μs  ┊ GC (min … max):  0.00% … 99.88%
 Time  (median):     505.769 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   583.345 ns ±   5.396 μs  ┊ GC (mean ± σ):  11.83% ±  4.49%

        ▃██▂                                                     
  ▁▁▁▂▃▆████▅▂▂▂▂▂▂▁▁▁▂▂▃▃▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  489 ns           Histogram: frequency by time          607 ns <

 Memory estimate: 880 bytes, allocs estimate: 6.

New code:

julia> @benchmark StrideArray(undef, real(mesh),
                              StaticInt(size(mesh.tree_node_coordinates, 1)),
                              static_length(solver.basis.nodes), static_length(mesh.nodes))
BenchmarkTools.Trial: 10000 samples with 142 evaluations.
 Range (min … max):  701.296 ns … 750.380 μs  ┊ GC (min … max):  0.00% … 99.87%
 Time  (median):     717.725 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   819.442 ns ±   7.505 μs  ┊ GC (mean ± σ):  10.81% ±  3.66%

    ▂▅▇█▇▆▅▂▁▁▂▃▂▃▃▄▄▄▃▃▃▂▂▂▂▂▁▁▁ ▁                             ▂
  ▆███████████████████████████████████▇▇▇▆▆▆▆▆▆▆▄▅▅▅▄▅▃▆▆▆▅▅▆▄▃ █
  701 ns        Histogram: log(frequency) by time        851 ns <

 Memory estimate: 880 bytes, allocs estimate: 6.

As we can see, here is some additional cost, but it seems to be a fixed overhead that doesn't depend on the mesh size, so I'm not sure if it's a significant regression in performance or not. I don't see any type instabilities from @code_warntype.

[ranocha]

Can you please also report benchmarks of a full call to calc_node_coordinates!?

[tristanmontoya]

Ok, actually there seems to be a type instability when I call the full calc_node_coordinates!, and the performance is worse. See below.

julia> using Trixi, BenchmarkTools

julia> using StaticArrayInterface: static_length

julia> using StrideArrays: StrideArray, StaticInt

julia> solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs);

julia> mesh = P4estMesh((8, 8), polydeg = 3,
                        coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0),
                        initial_refinement_level = 1);

julia> elements = Trixi.init_elements(mesh, LinearScalarAdvectionEquation2D((1.0,1.0)), solver.basis, Float64);

Benchmark using the existing code:

julia> @benchmark Trixi.calc_node_coordinates!(elements.node_coordinates, mesh, solver.basis.nodes)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  23.250 μs …  50.042 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     23.500 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   23.528 μs ± 487.328 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▁  ▆   █  ▆   ▅  ▄   ▁                                 
  ▂▁▁▃▁▁▁█▁▁█▁▁▁█▁▁█▁▁▁█▁▁█▁▁▁█▁▁▇▁▁▁▆▁▁▅▁▁▁▄▁▁▃▁▁▁▃▁▁▃▁▁▁▂▁▁▂ ▃
  23.2 μs         Histogram: frequency by time           24 μs <

 Memory estimate: 4.23 KiB, allocs estimate: 73.

Benchmark with the proposed change:

julia> @benchmark Trixi.calc_node_coordinates!(elements.node_coordinates, mesh, solver.basis.nodes)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  66.459 μs … 114.210 ms  ┊ GC (min … max):  0.00% … 99.92%
 Time  (median):     69.958 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   82.256 μs ±   1.142 ms  ┊ GC (mean ± σ):  14.26% ±  1.39%

    ▂▁      ▁     █▇▅ ▃▂▂                                       
  ▁▄██▇▃▂▂▄▇█▅▃▂▃▆███▇███▆▆▄▅▄▅▄▄▄▃▃▃▄▃▃▂▃▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  66.5 μs         Histogram: frequency by time         77.7 μs <

 Memory estimate: 45.66 KiB, allocs estimate: 784.

Because of these performance regressions, I will mark this PR as a draft until they are resolved, since we can work around this in TrixiAtmo.jl for now. Thanks @ranocha for the suggestions!

[ranocha]

By the way, you should create a new branch for PRs and not work from main of your fork

[tristanmontoya]

This will not be pursued further; instead, a new PR will be made with the ambient dimension stored as a type parameter.