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

High-order advection stencils, overhaul of boundary treatment, and bugfix in advection reconstruction #2603

Merged
merged 183 commits into from
Jul 25, 2022

Conversation

simone-silvestri
Copy link
Collaborator

@simone-silvestri simone-silvestri commented Jun 8, 2022

In this PR:

  • Refactoring of advection to have three different reconstruction schemes: Centered, UpwindBiased and WENO
  • advection schemes can be constructed with an order keyword argument which goes up to 12 for Centered and 11 for upwind schemes
  • Reimplementation of boundary treatment: reconstruction method _symmetric_, _left_biased_ and _right_biased have their own boundary treatment which entails checking the points for each method differently
  • advection scheme have a boundary_scheme field which will is used to reconstruct values near the boundary. Boundary schemes use the same method as the parent scheme but with a lower order:
    For example, when using a WENO 9th order, closing in near the boundaries the order will be progressively decreased to 7th, 5th, 3rd and finally 1st.

It also solves various bugs associated with advection and immersed boundary.

@navidcy navidcy added the numerics 🧮 So things don't blow up and boil the lobsters alive label Jun 8, 2022
@glwagner glwagner marked this pull request as draft June 9, 2022 14:38
@simone-silvestri
Copy link
Collaborator Author

I just realized that this might have been the problem which generated instabilities in #2510

@glwagner
Copy link
Member

glwagner commented Jun 10, 2022

I just realized that this might have been the problem which generated instabilities in #2510

Can you explain? How do you use the same interpolation schemes that we use for scalars to implement the term w * dz(u) ? I don't see instabilities on #2510 --- what were the results?

@simone-silvestri
Copy link
Collaborator Author

The problem was that we were generating spurious extrema when using high order advection near immersed boundaries. These spurious extrema were by no means negligible. I think that in the vertical direction this problem might have a larger impact on the stability due to the smaller grid scales and lower diffusion, while in the horizontal direction, diffusion could hide the extrema more effectively

const WENOVectorInvariantVel = Union{WENOVectorInvariantVel3, WENOVectorInvariantVel5}
const WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS, VI} =
Union{WENO3{FT, XT, YT, ZT, XS, YS, ZS, VI},
WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI}} where {FT, XT, YT, ZT, XS, YS, ZS, VI<:SmoothnessStencil}
Copy link
Member

@glwagner glwagner Jun 14, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is WENOVectorInvariant identified as a WENO with a SmoothnessStencil?

The design proposed by #2454 proposed to disentangle the reconstruction scheme from the "formulation":

struct VectorInvariantMomentumAdvection
     vorticity_stencil
     vorticity_reconstruction
end

Therefore "WENOVectorInvariant" is VectorInvariantMomentumAdvection with vorticity_reconstruction::WENO5.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn t get to that part yet

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay!

@glwagner
Copy link
Member

It doesn't look like things are moving towards #2454 exactly. Is there a plan to eventually do this, or is another design being proposed?

Copy link
Collaborator

@tomchor tomchor left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great work @simone-silvestri!

Tests seem to be passing so ready to merge?

@simone-silvestri
Copy link
Collaborator Author

simone-silvestri commented Jul 20, 2022

I d like to try to improve performance a bit first... Anyways, by next week I ll merge

@glwagner
Copy link
Member

@tomchor I guess we are hoping we can merge without the 3/4 performance penalty.

@simone-silvestri
Copy link
Collaborator Author

simone-silvestri commented Jul 22, 2022

comparison of run wall time for the bickley jet 128² between this PR and main (on the GPU).

This is running the code till 200 seconds

Periodic Periodic Bounded Bounded Immersed Immersed
Scheme this PR Main this PR Main this PR Main
Upwind 5 33.8967 64.0005 50.0858 47.4800 140.9554 107.9684
WENO 5 VI 28.4259 30.8202 37.8422 34.2834 123.7893 122.7207
WENO 5 FF 27.9952 31.9765 32.2018 36.7740 111.0842 107.0792
Vector Inv 28.9798 27.7792 37.5503 33.1336 125.4214 119.6501

Looking at these timings, everything seems in the ballpark (except for the immersed boundary) so I am not sure what might be the issue...

I guess we would need an in depth profiling to understand this?

I say we merge. The benefit from time step increase covers for the performance decrease. And finding the root of such a decrease might be very difficult

@glwagner
Copy link
Member

What's are the possible causes of the performance difference? What code should we focus on if we want to close this performance gap in the future?

@simone-silvestri
Copy link
Collaborator Author

simone-silvestri commented Jul 25, 2022

I believe it might be warp divergence. The way we have boundaries (and immersed boundaries) implemented, we create instances of branch divergence because the scheme is different near the boundaries. In this PR, for 5th order schemes there are two instances of divergence instead of one (the scheme changes twice instead of once). I am not super sure about it but it might be an explanation. A proper profile using a visual profiler will probably show it

@glwagner
Copy link
Member

Okay, so before this PR, we generate two stencils for each flux calculation --- second order and fifth order. The fifth order result is then used in the interior and the second order result near boundaries. In this new scheme, we generate three stencils rather than two stencils, because we limit to third order, and then to second order?

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>
@glwagner
Copy link
Member

For a future PR (probably #2642):

  • Use the order (not buffer size) as a type parameter, so WENO(order=5) yields a type WENO{5}, for example
  • Improve docstrings for advection schemes (but this is contingent on the larger advection refactor in Advection Refactor #2642)

@simone-silvestri simone-silvestri merged commit c2ac8fb into main Jul 25, 2022
@navidcy navidcy deleted the ss/advection_refactor branch June 18, 2023 22:46
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐞 Even a perfect program still has bugs numerics 🧮 So things don't blow up and boil the lobsters alive
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants