Add auxiliary variables

[benegee]

Add a container holding auxiliary node and surface node variables to semidiscretizations' cache, s.t. they can be used when computing rhs!.

As a first step, AcousticPerturbationEquations2DAuxVars are implemented. This is to demonstrate how auxiliary variables can be used instead of augmenting the solution state vector, as previously done, e.g., in AcousticPerturbationEquations2D.

This feature should ultimately be usable by

• [ ] @tristanmontoya: covariant solver on surfaces
• [x] @paulaweiss: variable coefficients for prescribed wind field
• [x] @benegee: perturbation approach for compressible Euler equations

To discuss:

• [x] naming: aux_vars vs. auxiliary_variables (and derived names such as get_auxiliary_surface_node_vars)
  The current plan is to use *aux*vars*

Closes #1623, closes #497, related #358

[github-actions[bot]]

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 :heart: by the Trixi.jl community.

[codecov[bot]]

Codecov Report

:x: Patch coverage is 54.91453% with 422 lines in your changes missing coverage. Please review. :white_check_mark: Project coverage is 95.80%. Comparing base (f4bbcd9) to head (19d0a72).

Files with missing lines	Patch %	Lines
src/solvers/dgsem_tree/dg_2d.jl	37.38%	201 Missing :warning:
src/solvers/dgsem_tree/containers_2d.jl	43.53%	96 Missing :warning:
src/equations/acoustic_perturbation_2d_aux_vars.jl	57.05%	64 Missing :warning:
src/solvers/dgsem_tree/dg_2d_parallel.jl	38.30%	29 Missing :warning:
src/visualization/utilities.jl	46.43%	15 Missing :warning:
src/solvers/dgsem_tree/indicators_2d.jl	46.15%	7 Missing :warning:
src/solvers/dgsem_tree/containers.jl	94.19%	5 Missing :warning:
src/solvers/dgsem_p4est/containers.jl	25.00%	3 Missing :warning:
src/callbacks_step/stepsize_dg2d.jl	90.00%	2 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2348      +/-   ##
==========================================
- Coverage   96.80%   95.80%   -1.00%     
==========================================
  Files         528      532       +4     
  Lines       42655    43524     +869     
==========================================
+ Hits        41292    41697     +405     
- Misses       1363     1827     +464     

Flag	Coverage Δ
unittests	95.80% <54.91%> (-1.00%)	:arrow_down:

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

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[sloede]

As discussed today, it's great to see this is coming to Trixi.jl! Feel free to ping me in case of questions or for a review!

[DanielDoehring]

@knstmrd is this maybe helpful for you?

[benegee]

CI now passes. Ready for comments @trixi-framework/developers

[benegee]

Maybe unrelated, but what was the idea behind https://github.com/trixi-framework/Trixi.jl/pull/576?

[sloede]

As discussed in the Trixi.jl meeting, it would be great if the first round of reviews could be done by @tristanmontoya 😊

[tristanmontoya]

As discussed in the Trixi.jl meeting, it would be great if the first round of reviews could be done by @tristanmontoya 😊

Yes, I will get to that soon!

[DanielDoehring]

This looks quite comprehensive - we should be really careful if this is breaking or not.

[jlchan]

Maybe unrelated, but what was the idea behind https://github.com/trixi-framework/Trixi.jl/pull/576?

I think I requested that, but I can't remember why now.

[benegee]

Finally, am I correct that the auxiliary variables are recomputed, rather than interpolated, when the mesh is adapted? Does this give a different result than adding them to the solution state u, since in the latter case they would be interpolated alongside the conservative variables?

This is correct. And in fact, this is certainly a crucial point. When, e.g., thinking about a hydrostatic background state, which is computed analytically, it makes sense to me to recompute it upon refinement. On the other hand, I can imagine that the induced numerical error in the system (interpolated solution and recomputed background state do not match anymore) could become relevant.

For the moment, I can only offer a simple comparison based on one example, which shows that both approaches agree.

With interpolated auxiliary variables (as part of the solution vector)

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'AcousticPerturbationEquations2D' with DGSEM(polydeg=5)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 88                run time:       3.98734360e-02 s
 Δt:             2.50000000e-01                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      3.00000000e+01 (100.000%)     time/DOF/rhs!:  2.82353752e-08 s
                                               PID:            3.28173783e-08 s
 #DOFs per field:          2736                alloc'd memory:        243.342 MiB
 #elements:                  76
 ├── level 5:                16
 ├── level 4:                36
 ├── level 3:                14
 └── level 2:                10

 Variable:       v1_prime         v2_prime         p_prime_scaled   v1_mean          v2_mean          c_mean           rho_mean
 L2 error:       1.94350488e-02   1.95225440e-02   4.81982254e-02   1.51483413e-16   0.00000000e+00   3.02966825e-16   3.02966825e-16
 Linf error:     1.83057645e-01   1.90845961e-01   1.03618809e+00   2.05391260e-15   0.00000000e+00   4.10782519e-15   4.10782519e-15
 ∑∂S/∂U ⋅ Uₜ :  -6.65307540e-08
────────────────────────────────────────────────────────────────────────────────────────────────────

With recomputed auxiliary variables:

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'AcousticPerturbationEquations2DAuxVars' with DGSEM(polydeg=5)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 88                run time:       2.44520200e-02 s
 Δt:             2.50000000e-01                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      3.00000000e+01 (100.000%)     time/DOF/rhs!:  1.70477983e-08 s
                                               PID:            2.00607297e-08 s
 #DOFs per field:          2736                alloc'd memory:        117.281 MiB
 #elements:                  76
 ├── level 5:                16
 ├── level 4:                36
 ├── level 3:                14
 └── level 2:                10

 Variable:       v1_prime         v2_prime         p_prime_scaled
 L2 error:       1.94350488e-02   1.95225440e-02   4.81982254e-02
 Linf error:     1.83057645e-01   1.90845961e-01   1.03618809e+00
 ∑∂S/∂U ⋅ Uₜ :  -6.65307540e-08
────────────────────────────────────────────────────────────────────────────────────────────────────

[amrueda]

Quoting @tristanmontoya

Finally, am I correct that the auxiliary variables are recomputed, rather than interpolated, when the mesh is adapted?

Quoting @benegee

This is correct. And in fact, this is certainly a crucial point.

I also think that this is an important point!! I can think of different reasons why you would like to have these variables interpolated/projected or recomputed. For instance:

• If you are dealing with a background state, maybe you do want to project the auxiliary variables if you are concerned about conservation. If you recompute them, the total mass/momentum/energy of the system might change after adapting the mesh.
• If you are dealing with exact geometric terms, as in Tristan's implementation of the covariant equations, you want the variables to be recomputed at the new points to ensure they correspond to the exact quantities.

Maybe in the future we would like to be able to select which variables get recomputed and which interpolated/projected. 🤔

[benegee]

Maybe in the future we would like to be able to select which variables get recomputed and which interpolated/projected. 🤔

Thanks for these examples @amrueda! It seems clear to me that we do want to be able to choose either approach.

[benegee]

I was also wondering, for problems without auxiliary variables, does this add any performance overhead? It doesn't look like it would be significant, but might be good to double check.

Yes, good point. So far, I did not see a significant difference in performance. Here is the output of elixir_euler_kelvin_helmholtz_instability.jl, just as an example.

main:

────────────────────────────────────────────────────────────────────────────────────
             Trixi.jl                      Time                    Allocations
                                  ───────────────────────   ────────────────────────
        Tot / % measured:              4.52s /  96.5%           50.9MiB /  99.4%

Section                   ncalls     time    %tot     avg     alloc    %tot      avg
────────────────────────────────────────────────────────────────────────────────────
rhs!                       3.15k    4.23s   97.1%  1.35ms   10.1KiB    0.0%    3.28B
  volume integral          3.15k    3.59s   82.5%  1.14ms      752B    0.0%    0.24B
    ~volume integral~      3.15k    3.32s   76.2%  1.06ms      752B    0.0%    0.24B
    blending factors       3.15k    272ms    6.2%  86.4μs     0.00B    0.0%    0.00B
  interface flux           3.15k    261ms    6.0%  83.1μs     0.00B    0.0%    0.00B
  prolong2interfaces       3.15k    183ms    4.2%  58.3μs     0.00B    0.0%    0.00B
  surface integral         3.15k    127ms    2.9%  40.4μs     0.00B    0.0%    0.00B
  reset ∂u/∂t              3.15k   33.6ms    0.8%  10.7μs     0.00B    0.0%    0.00B
  Jacobian                 3.15k   26.6ms    0.6%  8.45μs     0.00B    0.0%    0.00B
  ~rhs!~                   3.15k   5.24ms    0.1%  1.67μs   9.33KiB    0.0%    3.04B
  prolong2boundaries       3.15k    368μs    0.0%   117ns     0.00B    0.0%    0.00B
  prolong2mortars          3.15k    288μs    0.0%  91.4ns     0.00B    0.0%    0.00B
  mortar flux              3.15k    257μs    0.0%  81.6ns     0.00B    0.0%    0.00B
  boundary flux            3.15k    121μs    0.0%  38.4ns     0.00B    0.0%    0.00B
  source terms             3.15k   89.2μs    0.0%  28.3ns     0.00B    0.0%    0.00B
I/O                           34   66.6ms    1.5%  1.96ms   50.0MiB   98.9%  1.47MiB
  save solution               33   63.0ms    1.4%  1.91ms   49.9MiB   98.6%  1.51MiB
  get element variables       33   2.83ms    0.1%  85.7μs     0.00B    0.0%    0.00B
  ~I/O~                       34    838μs    0.0%  24.7μs    116KiB    0.2%  3.41KiB
  save mesh                   33   4.27μs    0.0%   129ns     0.00B    0.0%    0.00B
  get node variables          33   1.80μs    0.0%  54.6ns     0.00B    0.0%    0.00B
calculate dt                 630   36.5ms    0.8%  58.0μs     0.00B    0.0%    0.00B
analyze solution               8   23.4ms    0.5%  2.93ms    582KiB    1.1%  72.8KiB
────────────────────────────────────────────────────────────────────────────────────

This branch:

────────────────────────────────────────────────────────────────────────────────────
             Trixi.jl                      Time                    Allocations
                                  ───────────────────────   ────────────────────────
        Tot / % measured:              4.59s /  96.4%           50.9MiB /  99.4%

Section                   ncalls     time    %tot     avg     alloc    %tot      avg
────────────────────────────────────────────────────────────────────────────────────
rhs!                       3.15k    4.31s   97.4%  1.37ms   10.1KiB    0.0%    3.28B
  volume integral          3.15k    3.66s   82.8%  1.16ms      752B    0.0%    0.24B
    ~volume integral~      3.15k    3.38s   76.5%  1.08ms      752B    0.0%    0.24B
    blending factors       3.15k    276ms    6.2%  87.8μs     0.00B    0.0%    0.00B
  interface flux           3.15k    265ms    6.0%  84.2μs     0.00B    0.0%    0.00B
  prolong2interfaces       3.15k    186ms    4.2%  59.1μs     0.00B    0.0%    0.00B
  surface integral         3.15k    128ms    2.9%  40.7μs     0.00B    0.0%    0.00B
  reset ∂u/∂t              3.15k   34.5ms    0.8%  11.0μs     0.00B    0.0%    0.00B
  Jacobian                 3.15k   27.9ms    0.6%  8.86μs     0.00B    0.0%    0.00B
  ~rhs!~                   3.15k   4.44ms    0.1%  1.41μs   9.33KiB    0.0%    3.04B
  prolong2boundaries       3.15k    391μs    0.0%   124ns     0.00B    0.0%    0.00B
  prolong2mortars          3.15k    257μs    0.0%  81.7ns     0.00B    0.0%    0.00B
  mortar flux              3.15k    145μs    0.0%  46.1ns     0.00B    0.0%    0.00B
  boundary flux            3.15k    143μs    0.0%  45.5ns     0.00B    0.0%    0.00B
  source terms             3.15k   94.1μs    0.0%  29.9ns     0.00B    0.0%    0.00B
I/O                           34   54.5ms    1.2%  1.60ms   50.0MiB   98.9%  1.47MiB
  save solution               33   50.7ms    1.1%  1.54ms   49.9MiB   98.6%  1.51MiB
  get element variables       33   2.92ms    0.1%  88.4μs     0.00B    0.0%    0.00B
  ~I/O~                       34    938μs    0.0%  27.6μs    116KiB    0.2%  3.41KiB
  save mesh                   33   5.45μs    0.0%   165ns     0.00B    0.0%    0.00B
  get node variables          33   2.19μs    0.0%  66.4ns     0.00B    0.0%    0.00B
calculate dt                 630   36.8ms    0.8%  58.4μs     0.00B    0.0%    0.00B
analyze solution               8   23.5ms    0.5%  2.94ms    582KiB    1.1%  72.8KiB
────────────────────────────────────────────────────────────────────────────────────

[benegee]

Thanks for these examples @amrueda! It seems clear to me that we do want to be able to choose either approach.

I had a look.

• The aux vars vector in the cache is resized here
  https://github.com/trixi-framework/Trixi.jl/blob/4b7c029575fe027e0d43fefe209d138caee1c2d3/src/callbacks_step/amr_dg2d.jl#L119
• We would then need to call refine_element for the aux vars as well
  https://github.com/trixi-framework/Trixi.jl/blob/4b7c029575fe027e0d43fefe209d138caee1c2d3/src/callbacks_step/amr_dg2d.jl#L130
• Given the fact that something like if mesh isa P4estMesh is already used in the callback, it seems OK to me to use something like if aux_vars.aux_interpolation to distinguish between the interpolate and recompute cases.
• Similarly rebalance_solver! and coarsen! would need to be adapted.

Overall I think it is realistic to add this functionality. But I would postpone this to another PR.

[benegee]

I addressed #2386, #2387. If we postpone the interpolation approach to a subsequent PR, do you think this is ready for the next round of reviews? @amrueda @tristanmontoya

[tristanmontoya]

I addressed #2386, #2387. If we postpone the interpolation approach to a subsequent PR, do you think this is ready for the next round of reviews? @amrueda @tristanmontoya

Yes, I agree with postponing the interpolation approach and that the PR is ready for another round of reviews, but I think it should be documented somewhere (if it isn't already) that auxiliary variables get recomputed rather than interpolated.

[amrueda]

I addressed #2386, #2387. If we postpone the interpolation approach to a subsequent PR, do you think this is ready for the next round of reviews? @amrueda @tristanmontoya

Yes, I agree with postponing the interpolation approach and that the PR is ready for another round of reviews, but I think it should be documented somewhere (if it isn't already) that auxiliary variables get recomputed rather than interpolated.

Thanks @benegee! I also agree with postponing the interpolation and think that the PR is ready for the next round of reviews. 🙂

[benegee]

With #2400 and the latest changes in main merged, this is now ready for the next round of reviews.

Some comments:

• The auxiliary variables could potentially be used everywhere in Trixi.jl, which makes infrastructural changes necessary almost everywhere. It usually boils down to querying have_aux_node_vars(equations) and using the result for multiple dispatch.
• Specializations for have_aux_node_vars::True are only added to the 2D TreeMesh implementation in this PR (with others to follow). The code changes are trivial, usually just adding get_aux_node_vars or get_aux_surface_node_vars where get_node_vars or get_surface_node_vars were already called, and passing the resulting aux vars down, e.g. to flux computations.
• The only place where aux vars are actually used in the PR is equations/acoustic_perturbation_2d_aux_vars.jl, an alternative implementation of the already existing acoustic_perturbation_2d.jl equations.
• A less trivial part is the setup of the aux vars in solvers/dgsem_tree/containers_2d.jl, although it mostly follows corresponding prolong2* functions.
• A non-trivial part is the setup and usage of MPI interfaces, where the aux vars are directly evaluated on both sides and not communicated. This assume the aux vars field is steady and has no jumps.
• A truly non-trivial part is the setup and use of mortars. So far, aux vars are always evaluated for the small mortar elements and these values are used on both sides of the mortar, unlike the usual variables, which also exist on the large elements and need to be interpolated to the small elements. Similar to the discussion in https://github.com/trixi-framework/Trixi.jl/pull/2348#issuecomment-2862565759, an alternative interpolating version could be added.

@sloede @tristanmontoya @amrueda