Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIB: Subcell limiting with IDP and MCL approach #1502

Draft
wants to merge 558 commits into
base: main
Choose a base branch
from

Conversation

bennibolm
Copy link
Contributor

@bennibolm bennibolm commented Jun 1, 2023

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:

@bennibolm
Copy link
Contributor Author

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
Copy link
Member

sloede commented Jun 1, 2023

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
Copy link
Contributor Author

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
Copy link
Contributor Author

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
Copy link
Member

sloede commented Jun 1, 2023

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
Copy link
Contributor Author

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
Copy link
Member

sloede commented Jun 1, 2023

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
Copy link
Contributor Author

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
Copy link

codecov bot commented Jun 6, 2023

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
Copy link
Contributor Author

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
Copy link
Contributor Author

bennibolm commented Jul 6, 2023

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 and others added 18 commits February 5, 2024 13:56
* Use @Batch reduction for bounds check

* Remove section in documentation
* Merge main (first version)

* Fix bug

* Fix setup in tutorial

* Adapt tests

* Add bar_states variable to unit test

* fmt
* Add structured mesh support

* Fix non-periodic computation of bounds

* Use local limiting and nonperiodic domain in source terms elixir

* Use local limiting in free stream elixir

* Remove not needed lines

* Remove P4estMesh

* Add non-periodic tests with local bounds

* fmt

* Fix test

* Use `get_inverse_jacobian` instead of dispatching all routines

* Simplify `perform_idp_correction!`

* Revert stuff

* Remove free stream elixir

* Use sedov blast instead of source term setup; add news

* Update dispatching for mesh types

* Remove dispatching for mesh types

* Move new sedov tests within the test file

* Move new tests within test file

* Adapt dispatching

* Fix typo

* Remove not-needed parameter

* Add subcell limiting support for StructuredMesh (trixi-framework#1946)

* Add structured mesh support

* Fix non-periodic computation of bounds

* Use local limiting and nonperiodic domain in source terms elixir

* Use local limiting in free stream elixir

* Remove not needed lines

* Remove P4estMesh

* Add non-periodic tests with local bounds

* fmt

* Fix test

* Use `get_inverse_jacobian` instead of dispatching all routines

* Simplify `perform_idp_correction!`

* Revert stuff

* Remove free stream elixir

* Use sedov blast instead of source term setup; add news

* Update dispatching for mesh types

* Move new tests within test file

* Adapt dispatching

* Fix typo

* Remove not-needed parameters

* Dispatch `check_bounds` for dimension using u
@bennibolm
Copy link
Contributor Author

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.

bennibolm and others added 4 commits September 20, 2024 13:52
* Add tests for alpha analysis

* Add comments; Fix two tests

* Fix test

* Update alpha tests

* Rm not needed code

* Simplify test
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants