WIB: Subcell limiting with IDP and MCL approach

[bennibolm]

This PR includes the complete implementation of subcell limiting methods (IDP and MCL). The volume integral VolumeIntegralSubcellLimiting is used. This contains the following features:

IDP Limiting: SubcellLimiterIDP (uses stage callback SubcellLimiterIDPCorrection)

• Local two-sided Zalesak-type limiter for conservative variables local_twosided_variables_cons
• Positivity limiter for positivity_variables_cons and positivity_variables_nonlinear
• One-sided limiter for non-linear variables local_onesided_variables_nonlinear

MCL: SubcellLimiterMCL

• density_limiter: Local maximum/minimum limiter for cons(1)
• density_coefficient_for_all: Apply cons(1) blending coefficient for all quantities
• sequential_limiter: Local maximum/minimum for variables phi:=cons(i)/cons(1) for i in 2:nvariables
• conservative_limiter: Local maximum/minimum for conservative variables 2:nvariables
• positivity_limiter_pressure: Positivity limiter for pressure à la Kuzmin. (Sharp version with positivity_limiter_pressure_exact=true
• positivity_limiter_density: Positivity limiter for cons(1)
• entropy_limiter_semidiscrete: Semi-discrete entropy fix

Additional features:

• Bounds Check routine in every RK stage (stage callback)
• Export of blending coefficients (node variables) in the SaveSolutionCallback
• Analysis of blending coefficient for every time step (step callback)
• Smoothness indicator IndicatorHennemannGassner (hard switch)

TODOs:

• Merge stage callback use of antidiffusive stage into Subcell positivity IDP limiting for conservative variables #1476, if desired
• Merge Multicomponent Euler elixir and test into Subcell positivity IDP limiting for conservative variables #1476
• Add tests using StructuredMesh
• Generalize density_tvd for arbitrary conservative variables (Generalizing two-sided state limiting for conservative variables bennibolm/Trixi.jl#112)
• Relocate resize!(semi, nelements) to allow using Time Integrations methods by OrdinaryDiffEq for MCL(?):
  - Possible at end of semidisretize?
  - in Elixirs before solve function?
• Fix export of node variables for parallel meshes (causing failing MPI tests) (still a problem?)
• Update names of parameters for MCL using snake_case
• Merge PR of new analysis routine of blending factors (Implement alpha analysis as Callback bennibolm/Trixi.jl#110)
• Implement goal functions for Newton method for arbitrary variables and equations (Support nonlinear positivity limiting for arbitrary variables bennibolm/Trixi.jl#115)
• Move auxiliary IDP and MCL routines to new file dg_2d_subcell_limiters.jl, subcell_limiters(_2d).jl? Issue:
  • dg_2d and indicators_2d file much longer than before. -> Created specific files for subcell limiter code
• Implement variable_bounds as Dict? Saves a lot of messy stuff. Subcell limiting: Revise order of bounds using a Dict #1649
• Implement MCL entropy fix for non-Cartesian meshes
• Reduce simulation time of tests
• Add comment explaining why we implemented the SSPRK scheme with numerator and denominator. Link PR Assure conservation for SSP scheme #1640
• Possible to add Limiting for coupled semidiscretizations and other semidiscretizations?
• Nonconservative support
• Is create_cache() type stable when asking if bar_sates ... end?
• Plotting variable: Keeping, deleting, renaming?
• Thread parallelize bounds check for IDP (Thread-)Parallelize bounds check routine for subcell IDP limiting #1736
• Thread parallelize bounds check for MCL Parallel bounds check for MCL bennibolm/Trixi.jl#124
• IDP Limiting with P4estMesh Add support for P4estMesh to subcell-limiting bennibolm/Trixi.jl#121
• General one-sided local limiting for non-linear variables Generalize one-sided limiting bennibolm/Trixi.jl#125
• Clean up branch (delete unnecessary elixirs, clean up tests)

[bennibolm]

This is currently the newest and most complete version of Subcell Limiting. The different limiters are briefly mentioned above.

In addition, there are some smaller TODOs, but they don't have much to do with the first (already existing) IDP positivity limiter #1476. For this one, I would only add the implementation of the antidiffusive stage as a stage callback. This allows other stage callbacks in the future (as the BoundsCheck routine or possible IMEX-type-Methods for Multi-Ion simulations, @amrueda is thinking about).
All other new features can be added bit by bit.

So the question. Do we keep the current structure with its own volume integral VolumeIntegralShockCapturingSubcell (perhaps with a different name since it is no longer just shock capturing) or do we change this. @sloede @gregorgassner

[sloede]

So the question. Do we keep the current structure with its own volume integral VolumeIntegralShockCapturingSubcell (perhaps with a different name since it is no longer just shock capturing) or do we change this. @sloede @gregorgassner

Can you elaborate a little bit? That is, what are the alternatives and what are the pros/cons for the various approaches?

[bennibolm]

Can you elaborate a little bit? That is, what are the alternatives and what are the pros/cons for the various approaches?

I personally think that the current structure is not bad. Rn we need some dispatches to distinguish between VolumeIntegralShockCapturingSubcell and other volume integrals. But this is actually the only “disadvantage”.

An alternative would be a custom solver for the subcell limiting. I think this strict distinction would lead to more confusion than to advantages (if there are any at all). Especially after the changes from the last week, where one is very flexible and can also run Elixirs without limiting with the implemented SSPRK algorithm.

[bennibolm]

I personally think that the current structure is not bad. Rn we need some dispatches to distinguish between VolumeIntegralShockCapturingSubcell and other volume integrals. But this is actually the only “disadvantage”.

An alternative would be a custom solver for the subcell limiting. I think this strict distinction would lead to more confusion than to advantages (if there are any at all). Especially after the changes from the last week, where one is very flexible and can also run Elixirs without limiting with the implemented SSPRK algorithm.

I have to revise myself a bit. Due to the fact that we partly use limiters based on the so-called bar states, we need the boundary condition and the current time within the volume integral for their calculation.
So, rn we dispatch calc_volume_integral! (see here). This works, but does not really correspond to the separation of the different integral functions, as it is currently valid in Trixi.

This could not really be solved by a separate solver, but the problem would then only exist there.

[sloede]

So, rn we dispatch calc_volume_integral! (see here). This works, but does not really correspond to the separation of the different integral functions, as it is currently valid in Trixi.

I don't see this additional dispatch in the current PR? Is there a third PR coming? Also, it would be great to clarify the relationship between #1476 and this one here.

[bennibolm]

I don't see this additional dispatch in the current PR? Is there a third PR coming? Also, it would be great to clarify the relationship between #1476 and this one here.

As I said, some limiter uses the bar states to calculate the bounds. The positivity limiting of IDP limiting (as included in the current PR) does not require the bar states. That's why I omitted it for now.

I thought it would be the best to divide this big PR into many small ones and started with the positivity limiter. When it is merged, there are some smaller PRs with simple additional routines. In total, there will be many more to come. I will add an overview over the PR's relationships in the description. Hopefully this makes it clearer.

[sloede]

I thought it would be the best to divide this big PR into many small ones and started with the positivity limiter. When it is merged, there are some smaller PRs with simple additional routines. In total, there will be many more to come. I will add an overview over the PR's relationships in the description. Hopefully this makes it clearer.

Ah, never mind, I had looked at this in the wrong way, and I now found the relevant sections where the calc_volume_integral! routine is changed... Phew. I think this is is something to be discussed with the others, maybe after the next Trixi Meeting?

From a first look, the dg_2d.jl file gets more than doubled in number of lines by this addition. That could motivate to make it a whole new solver "family", since it is becoming very far removed from a classical DG method. OTOH, since you really only modify the volume integral, it fits were well into the existing solver structure and creating a new solver for it might create more headaches than it helps - it could be viewed as a similar extension to the purely hyperbolic DG scheme as the parabolic stuff is.

Thus, for now I tend to lean towards keeping it as a super-special volume integral but to move everything IDP/MCL-related to a new file dg_2d_idp_mcl.jl (or similar) and to include that at the bottom of src/solvers/dgsem_tree/dg.jl

[bennibolm]

Phew. I think this is is something to be discussed with the others, maybe after the next Trixi Meeting?

Yeah, I know, that's heavy. Maybe it would be better to remove it from the volume integral. But then, we would have an additional function somewhere.

From a first look, the dg_2d.jl file gets more than doubled in number of lines by this addition. That could motivate to make it a whole new solver "family", since it is becoming very far removed from a classical DG method. OTOH, since you really only modify the volume integral, it fits were well into the existing solver structure and creating a new solver for it might create more headaches than it helps - it could be viewed as a similar extension to the purely hyperbolic DG scheme as the parabolic stuff is.

Thus, for now I tend to lean towards keeping it as a super-special volume integral but to move everything IDP/MCL-related to a new file dg_2d_idp_mcl.jl (or similar) and to include that at the bottom of src/solvers/dgsem_tree/dg.jl

I wasn't sure if it all made sense in the current location anyway. So, I think it's a good idea.

[codecov]

Codecov Report

Attention: Patch coverage is 95.08434% with 102 lines in your changes missing coverage. Please review.

Project coverage is 96.30%. Comparing base (4a81e00) to head (c7263f7).

Files with missing lines	Patch %	Lines
src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl	92.30%	58 Missing ⚠️
src/callbacks_step/stepsize_dg2d.jl	67.24%	19 Missing ⚠️
src/callbacks_step/limiting_analysis.jl	92.86%	6 Missing ⚠️
src/callbacks_step/limiting_analysis_2d.jl	93.42%	5 Missing ⚠️
src/solvers/dgsem_tree/subcell_limiters.jl	95.18%	4 Missing ⚠️
src/equations/compressible_euler_2d.jl	94.55%	3 Missing ⚠️
src/time_integration/methods_SSP.jl	91.18%	3 Missing ⚠️
...solvers/dgsem_structured/dg_2d_subcell_limiters.jl	98.57%	2 Missing ⚠️
...es/structured_2d_dgsem/elixir_euler_double_mach.jl	96.67%	1 Missing ⚠️
...tructured_2d_dgsem/elixir_euler_double_mach_MCL.jl	96.67%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1502      +/-   ##
==========================================
- Coverage   96.36%   96.30%   -0.06%     
==========================================
  Files         477      500      +23     
  Lines       37753    39791    +2038     
==========================================
+ Hits        36378    38319    +1941     
- Misses       1375     1472      +97     

Flag	Coverage Δ
unittests	96.30% <95.08%> (-0.06%)	⬇️

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.

[bennibolm]

So, finally we had our discussion yesterday about the naming issue. To summarize the results, I think we could agree on the names SubcellLimiterIDP and SubcellLimiterMCL for the current indicators. Then, I define these limiters by limiter_idp and limiter_mcl to include them to VolumeIntegralSubcellLimiting.
We also had a few suggestions for the a posteriori antidiffusive stage of SubcellLimiterIDP (implemented as a stage callback). I personally like for instance SubcellLimiterIDPCorrection to make the connection to SubcellLimiterIDP clear.

Additionally, I will add all possible remarks (in the documentation and the code itself) for IDP limiting, to clarify that SubcellLimiterIDP and SubcellLimiterIDPCorrection each do not work alone.

Is everyone happy with this or does anyone have any other suggestions? @sloede @ranocha @amrueda @gregorgassner @andrewwinters5000 @jlchan

[bennibolm]

I was thinking about the reordering of the subcell limiting code. First, of course, the new limiters can't stay in indicators.jl. I moved them to a new file named subcell_limiters.jl. Accordingly, there are now also subcell_limiters_2d.jl.
Then, moreover, there were many routines located in dg_2d.jl. I moved them to a new dg file dg_2d_subcell_limiters.jl.
This restructuring also applies for the respective files for StructuredMesh.
The definition of all containers (there are 4 in total) are still in containers_2d.jl.

An alternative would be an additional folder dgsem_idp_mcl (or dgsem_subcell_limiting). This would have the disadvantage that the implementation for all meshes would be together in there. (Or we add one folder for each mesh, which would certainly be too much.) Another option: We just use one new folder, and only the files distinguish between the meshes, e.g., dg_tree_2d.jl. However, this would be completely contrary to the existing structure of Trixi.

I prefer the current option. What do you think @sloede @ranocha? Maybe it would be useful to discuss this as well. Possibly in a smaller round.

[bennibolm]

In this discussion, I thought that I could easily remove mesh as a parameter of get_boundary_outer_state. Now, after merging those changes into this big subcell-limiting branch, I realize that it is not that easy as I thought. Because...

When such a function is called within a simulation using a StructuredMesh the parameters (see here) are

u_inner, t, boundary_conditions[1], Ja1, 1, equations, dg, cache, 1, j, element

Since the function header of get_boundary_outer_state includes indices... the number of parameters is not fixed.
So, when having multiple versions of get_boundary_outer_state for different mesh types as wanted in the future, there could be an error due to multiple function candidates or a completely wrong choice.
Because of that, for now, I kept the parameter mesh at least in this branch here.

There is already an issue about the interface of get_boundary_outer_state (see #2053). So, I will not add another, but wanted to formulate the problem somewhere.