Add capability to handle with node variables

[bennibolm]

This PR adds the functionality to handle node-level data beside the solution itself.
This is required to visualize subcell limiting coefficients, as used in #1476 and #1502.

It comes together with the PR to Trixi, which saves the coefficients in the .h5 files.

TODO:

• Adapt test to "new" testing structure (compare value and not only check whether the code runs through). Check this after the Trixi PR is merged.

Question:

• Useful to allow reinterpolation for subcell blending coefficients? Actually not, but it would be very annoying to always handle two files. Added a comment.

[bennibolm]

The current implementation supports both reinterpolation and no reinterpolation for the subcell limiting coefficient. That's the way it is handled for the solution itself.
Does it even make sense to allow reinterpolation for this coefficient? What do you think? @amrueda @andrewwinters5000 @sloede

I already prepared an PR to Trixi2Vtk_reference_files to add the needed .vtu files for testing. I will open the PR when this discussion here has been decided.

[sloede]

Does it even make sense to allow reinterpolation for this coefficient?

IMHO, no. It's like an AMR marker (should this cell be refined or not), which is not a continuous field that happens to be represented by a piecewise-constant approximation. Thus reinterpoaltion does not give a meaningful representation of reality but can actually be harmful.

Having said that, I am not sure how easy/annoying it is to always create two files, once once with reinterpolation (for the solution itself) and once without (for the node variables). It might thus be nice to have the ability for reinterpolation, just so you can have only one VTK file for a quick'n'dirty visualization.

What we really need/want though, is a ParaView plugin that can read in our HDF5 natively without us having to convert to VTK first. But that's for another day I guess :-)

[amrueda]

It might thus be nice to have the ability for reinterpolation, just so you can have only one VTK file for a quick'n'dirty visualization.

I agree with @sloede

[bennibolm]

Does it even make sense to allow reinterpolation for this coefficient?

IMHO, no. It's like an AMR marker (should this cell be refined or not), which is not a continuous field that happens to be represented by a piecewise-constant approximation. Thus reinterpoaltion does not give a meaningful representation of reality but can actually be harmful.

I thought so too. But I implemented the reinterpolation of the coefficients before it was possible to not reinterpolate the solution and just kept the code after the update.

Having said that, I am not sure how easy/annoying it is to always create two files, once once with reinterpolation (for the solution itself) and once without (for the node variables). It might thus be nice to have the ability for reinterpolation, just so you can have only one VTK file for a quick'n'dirty visualization.

I'm not sure that I understand your thoughts correctly. I don't think we need two different files here either way. It should be enough to just delete the implication of parameter reinterpolation for the subcell stuff.

EDIT: Oh, I now see the problem. I didn't realized that it's a problem to save reinterpolated and non-reinterpolated node-level data in one .vtu file.

What we really need/want though, is a ParaView plugin that can read in our HDF5 natively without us having to convert to VTK first. But that's for another day I guess :-)

😌

[amrueda]

I'm not sure that I understand your thoughts correctly. I don't think we need two different files here either way. It should be enough to just delete the implication of parameter reinterpolation for the subcell stuff.

No, I don't think that's a good idea. In some situations, it might still be advantageous to reinterpolate the solution variables to get a nice and smooth representation of the (part of the field that is a pure) DG solution.

[bennibolm]

Having said that, I am not sure how easy/annoying it is to always create two files, once once with reinterpolation (for the solution itself) and once without (for the node variables). It might thus be nice to have the ability for reinterpolation, just so you can have only one VTK file for a quick'n'dirty visualization.

I'm not sure that I understand your thoughts correctly. I don't think we need two different files here either way. It should be enough to just delete the implication of parameter reinterpolation for the subcell stuff.

No, I don't think that's a good idea. In some situations, it might still be advantageous to reinterpolate the solution variables to get a nice and smooth representation of the (part of the field that is a pure) DG solution.

After understanding the problem correctly 😅 , I will keep the current implementation for now and add a warning in that case, if that's okay for you.

[andrewwinters5000]

I also think that it not very "meaningful" to reinterpolate the node-wise alpha values because it could lead to misinterpretation (as already pointed out by the others). For the element-wise variant we can get away with this re-interpolation because the alphas are constant within each element, so it does not change the values.

Having said that, I am not sure how easy/annoying it is to always create two files, once once with reinterpolation (for the solution itself) and once without (for the node variables). It might thus be nice to have the ability for reinterpolation, just so you can have only one VTK file for a quick'n'dirty visualization.

I'm not sure that I understand your thoughts correctly. I don't think we need two different files here either way. It should be enough to just delete the implication of parameter reinterpolation for the subcell stuff.

I think that @sloede is correct. It would be annoying but you could create two separate files where the solution re-interpolates but the alphas do not. It would be something like (note I have not tested this)

trixi2vtk(joinpath("out", "solution_*.h5"), output_directory="out")
trixi2vtk(joinpath("out", "solution_*_celldata.h5"), output_directory="out", reinterpolate=false)

The second call might need to use ? instead of the * for the different files numbers.

[bennibolm]

So, to avoid the annoying part with two files, I just added a comment if reinterpolation=true is called with a subcell coefficient.

For the testing, I added two new tests. After running elixir_euler_shockcapturing_subcell.jl I check for reinterpolation=true and false. I would suggest, that I first request a review on the current changes here and then, after every is okay here, I open a PR to Trixi2Vtk_reference_files to add the files and (hopefully) fix both new tests.

[sloede]

Do I get this correctly that we don't run the "normal" tests on macos and windows? Only the threaded and mpi tests.

Yes. Since fewer macOS runners are available to us on the free plan, we only run some basic tests on macOS and assume that everything else "just works" as well. Originally, we tested everything everywhere, but this just wasn't feasible anymore due to long CI runtimes.

Do you think I sould spend some time on investigating this? I mean, the method itself is very sensitive to small changes by construction.

This is something you need to discuss with your advisor, if it is worth to investigate right now (maybe bring it up on Friday?). IMHO it is at least somewhat unfortunate that floating point truncation errors seem to have such a significant impact already. It indicates to me that either there's something wrong with the method itself (unlikely), the implementation (less unlikely, but still), or that the setups are not robust enough (my current favorite). By robust I mean that maybe you are not running a stable simulation but rather simulations that are just not crashing. I wonder if choosing slightly different parameters would allow you to alleviate the issue (thus confirming my suspicion) or whether it is really by design.

A short-term fix for this current PR here, could be to accept the differences and exclude the test from ci for macOS.

Yes. But please leave a short comment explaining the reason for this and x-ref this PR (#67)

[bennibolm]

IMHO it is at least somewhat unfortunate that floating point truncation errors seem to have such a significant impact already. It indicates to me that either there's something wrong with the method itself (unlikely), the implementation (less unlikely, but still), or that the setups are not robust enough (my current favorite). By robust I mean that maybe you are not running a stable simulation but rather simulations that are just not crashing. I wonder if choosing slightly different parameters would allow you to alleviate the issue (thus confirming my suspicion) or whether it is really by design.

Yes, I agree. I took a quick look into it and made a list of failing tests. I will add this here but will present this also on friday:
Subcell Limiting in Trixi:

• Failing tests for julia 1.11
    ◦ elixir_euler_blast_wave_sc_subcell_nonperiodic.jl 
      initial_condition_blast_wave
      local limiting of rho and entropy_math
      initial refinement level = 4 (in test)
      Dirichlet Boundary Conditions
    ◦ elixir_euler_sedov_blast_wave_sc_subcell.jl
      initial_condition_sedov_blast_wave
      local limiting of rho and entropy_guermond_etal
      initial_refinement_level=4
      periodic
    ◦ elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl
      initial_condition_shock_bubble
      local limiting of rhos
      initial refinement level = 3
      periodic
• Working tests:
    ◦ elixir_euler_shockcapturing_subcell.jl (Positivity of rho, initial_refinement_level = 5)
    ◦ elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl (Positivity of rho and p, in_refine_level=5)
    ◦ elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity (Positivity of rhos, in_ref_level=3)
    ◦ elixir_mhd_shockcapturing_subcell (Positivity of rho and p, initial_refinement_level=4)
    ◦ elixir_euler_sedov_blast_wave_sc_subcell (StructuredMesh, Positivity of Rho and p, 16x16)
    ◦ elixir_euler_sedov_blast_wave_sc_subcell (StructuredMesh, Local limiting of rho and entropy_guermond, 16x16)
• Used setup in tests of Trixi2Vtk
    ◦ elixir_euler_sedov_blast_wave_sc_subcell.jl (see above, local limiting, within failing tests for switch to julia 1.10 and 1.11)
• Switch to julia 1.10
    ◦ elixir_euler_sedov_blast_wave_sc_subcell.jl
    ◦ elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl

After that (but of course this must be investigated further) I would assume your are right with:

[...] or that the setups are not robust enough (my current favorite). By robust I mean that maybe you are not running a stable simulation but rather simulations that are just not crashing. I wonder if choosing slightly different parameters would allow you to alleviate the issue (thus confirming my suspicion) [...]

[bennibolm]

After stabilizing some elixirs in trixi-framework/Trixi.jl#2007 and a new Trixi release, I'm about to update the reference files in the specific repository (trixi-framework/Trixi2Vtk_reference_files#4) and then, hopefully, we can merge this PR 😅

[sloede]

LGTM!

[sloede]

@bennibolm do you know why the tests do not seem to run?

[bennibolm]

@bennibolm do you know why the tests do not seem to run?

@sloede I think, they are running, but are somehow showed very strange. When I press on the red X at the last commit, the CI test ran and there is an old version (failing) and a new version of the coverall tests.

I will try to make them run and ping you again.

[bennibolm]

@sloede I adapted the test setup for macOS. Is it possible that one needs to adapt the tests in GitHub that must pass in a PR? That could be the reason why they are not shown correctly but actually run properly.

[sloede]

I don't understand why the tests don't show up here, but as you pointed out correctly, they seem to have passed for the latest commit. Thus I'll merge this now.

[bennibolm]

Thanks