MPData advection scheme #3434
MPData advection scheme #3434
Conversation
simone-silvestri
added
🚨 DO NOT MERGE 🚨
| architecture :: Arch | ||
| underlying_grid :: G | ||
| immersed_boundary :: I | ||
| interior_active_cells :: M | ||
|
|
||
| surface_active_cells :: S |
There was a problem hiding this comment.
what surface is this, top, bottom, north, south?
There was a problem hiding this comment.
That is meant to be the top boundary. It comes from PR #3404. I have to correct it there.
There was a problem hiding this comment.
Actually, to be precise, it's the active map of reduced fields in the z-direction
There was a problem hiding this comment.
That is only relevant for the hydrostatic model right? It wouldn't apply to the non-hydrostatic model, or probably other models that are isotropic and 3D.
|
@simone-silvestri Have you already done tests to ensure that it is positive preserving? |
|
Indeed, it is positivity preserving. Not bounds preserving though. And it seems to be too diffusive. |
That's a big step forward! I know there are ways to iterate to improve the accuracy. I presume your test was the low order option? Do you think we should try more iterations to improve the upper bound and general accuracy? |
|
I have introduced the option to iterate as well an option to use the formulation that assumes infinite iterations ( example: grid = RectilinearGrid(size = (10, 10, 10), extent = (1, 1, 1))
advection = MPData(grid; iterations = 1) # equivalent to UpwindBiased(; order = 1)
advection = MPData(grid; iterations = 3) # 3 corrective iterations
advection = MPData(grid; iterations = nothing) # ∞ iterations (a little more expensive) |
|
This is the output of the 1D advection test in
|
|
with a 10 times smaller time step |
| @@ -104,6 +105,8 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac | |||
| calculate_pressure_correction!(model, first_stage_Δt) | |||
| pressure_correct_velocities!(model, first_stage_Δt) | |||
|
|
|||
| correct_advection!(model, Δt) | |||
There was a problem hiding this comment.
that it is only upwind first order, definitely incorrect 😄
There was a problem hiding this comment.
maybe "refine advective tendencies" is the right word?
|
Thanks for sharing the plots @simone-silvestri! There does not seem to be much improvement between 2 and 5 iterations. In the coarser case they look the same and in the finer case the more iterations has a higher maximum. The method is not supposed to preserve a maximum is it? |
|
Nono, it is not bounds preserving. It is only intended for positive definite quantities and prevents negative values. |
|
I produced a few more figures that might be of interest. First, it's the same plot as @simone-silvestri produced but with fifth order upwinding.
|
|
The second is with double the spatial resolution, and therefore halving the time step. |
|
The third and final one is with doubling the spatial resolution again (total of x4) and calculating the time step in the same way as usual.
|
|
Some thoughts: The first case ( The second case ( The third case ( I lesson that I learned from here is to use as much spatial resolution as we can as this does much more compared to the number of iterations. @simone-silvestri : do you think we should do a similar comparison for 2D advection? |
|
Thanks for running these tests. It would be nice to see how MPData performs in a 2D flow. It's interesting that tweaking the number of iterations doesn't seem to make much difference. I'm a bit skeptical about this scheme though; it seems pretty diffusive and doesn't even beat third-order upwind which is way simpler. I might have introduced a bug although I went through the code again and it looks correct to me. Another possibility is that the high diffusivity is the trade-off for keeping the scheme positivity-preserving? |
|
I looked into MPData a long time ago and I recall it being very diffusive, so it might be the nature of the beast. It beats first order upwinding, which is also positive preserving, but not by much. The WENO schemes seem to be a more complicated way to ensure positivity more accurately but I remember you saying there were problems doing this properly in multiple dimensions. There is no perfect scheme for sure. This is a paper that presents a positive preserving WENO scheme for population models. Example 3 is a 2D problem but the plots are only 1D. How similar is this to what you have tried already? |
|
This paper presents a maximum principle preserving flux limiter for FD RK-WENO schemes and has applications to incompressible flows in 2D. They claim to maintain accuracy and not have as small a time step as other approaches. I'm not sure of the details thought but I thought I would share, in case it was of interest. |
This PR implements a three-dimensional positivity preserving advection scheme called multidimensional positive definite advection transport algorithm (MPDATA), which follows this paper.
We don't need to merge this PR, but I'll leave it here in case, in the future, we need a multidimensional positivity-preserving advection scheme (for example for wetting and drying).