Maximum-principle-satisfying WENO scheme #2557

simone-silvestri · 2022-05-16T21:57:39Z

Following https://royalsocietypublishing.org/doi/epdf/10.1098/rspa.2011.0153 (section 5)

At the moment this is just a scheme check, so I have implemented it in 1D (and only in the horizontal direction)

For the moment this script

using Oceananigans

arch = CPU()
grid = RectilinearGrid(arch, size = (100, 100, 1), x = (0, 2π), y = (0, 2π), z=(0, 1), topology = (Periodic, Periodic, Bounded))

free_surface = ImplicitFreeSurface()

model = HydrostaticFreeSurfaceModel(; grid, free_surface,
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(bounds = (0, 1)),
                                    buoyancy = nothing,
                                    tracers = :c)

@inline function linear_u_profile(x, y, z) 
    if x < π/2 || x > 3 * π/2
        return 0.0
    else
        if y > 7/4*π || y < π/4
            return 0.0
        elseif y < 5/4 * π && y > 3/4 * π
            return π/2
        elseif y >= π/4 && y <= 3/4 * π
            return y - π/4
        else
            return π/2 - (y - 5/4 * π)
        end
    end
end

@inline function rectangle_c(x, y, z)
    if x > 2/3 * π && x < 4/3 * π
        if x > 2/3 * π && x < 4/3 * π
            return 1.0
        else
            return 0.0
        end
    else
        return 0.0
    end
end

set!(model, u = linear_u_profile, c = rectangle_c)

Δt = grid.Lx / grid.Nx / maximum(abs, model.velocities.u) * 0.2

for step in 1:10000
    time_step!(model, Δt)
    @info "step $step, time $(model.clock.time)"
end

returns

julia> model.tracers.c
100×100×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 100×100×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 106×106×7 OffsetArray(::Array{Float64, 3}, -2:103, -2:103, -2:4) with eltype Float64 with indices -2:103×-2:103×-2:4
    └── max=1.03919, min=-0.0214752, mean=0.351983

for usual WENO5
and

julia> model.tracers.c
100×100×1 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 100×100×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 106×106×7 OffsetArray(::Array{Float64, 3}, -2:103, -2:103, -2:4) with eltype Float64 with indices -2:103×-2:103×-2:4
    └── max=0.999906, min=1.09505e-5, mean=0.352073

for WENO5(bounds = (0, 1))

So it seems like this scheme, even in its 1D formulation, might help with our tracer bound issues.

simone-silvestri · 2022-05-16T22:28:56Z

This is a heatmap of the tracer after 10000 time steps (the initial condition of the tracer is 0 <= c <= 1). The top is WENO5, the bottom WENO5(bounds=(0, 1)). The out-of-bounds cells are substituted with a NaN

simone-silvestri · 2022-05-16T23:45:40Z

Same plot obtained with a second order centered scheme for tracers

francispoulin · 2022-05-17T00:43:24Z

This is great to see!

I imagine we could use this for the height field in ShallowWaterModel to allow for wetting and drying for a single layer. And also to prevent layer collapse in the case of multiple layers.

simone-silvestri · 2022-05-17T00:50:03Z

yeah, I guess in that case you don't want to set an upper bound so we can use WENO5(bounds = (0, Inf))

glwagner · 2022-05-17T02:50:57Z

Just a small plea to make this implementation compatible with / moves us towards the refactor envisioned in #2454

francispoulin · 2022-05-17T12:57:33Z

src/Advection/positivity_preserving_tracer_advection_operators.jl

+    return 1/Vᶜᶜᶜ(i, j, k, grid) * (div_x + div_y + div_z)
+end
+
+@inline function bounded_tracer_flux_divergence_x(i, j, k, grid, advection::BoundPreservingScheme, u, c)


This looks great! It seems like these three funcionts are almost identical except for u -> v -> w. Is that true?

If yes, I wonder if there is away to merge these to one, rather than having three functions that do the same thing?

Indices, areas and interpolations differ. I can try merging it with some metaprogramming

Not necessary. This is maybe cleaner but thought I would make the observation.

I prefer having separate functions over metaprogramming in this case.

@francispoulin the difference between x, y, z is that they depend on indices i, j, k respectively. So we also have advective_tracer_flux_x, etc.

glwagner · 2022-05-17T14:48:25Z

src/Advection/weno_fifth_order.jl

@@ -26,7 +30,7 @@ Weighted Essentially Non-Oscillatory (WENO) fifth-order advection scheme.

 $(TYPEDFIELDS)
 """
-struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} <: AbstractUpwindBiasedAdvectionScheme{2}
+struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF, PP} <: AbstractUpwindBiasedAdvectionScheme{2}


Is this the design we will want when we do #2454 ? If we can avoid digging ourselves into a deeper hole we should, but maybe it's not possible without a lot of effort...

I am not still completely sure how #2454 will go along, this change only modifies the div_Uc so I guess it might be easy to adapt in the future

In #2454 we will disentangle vector invariant from WENO5, so VI will go away. It might still make sense to blend WENO5 with flux-limited WENO5, however. It's just that if we have nth order WENO, we would also have to figure out how to do nth order flux limited WENO if they are entangled.

I think the limiter will be independent on the order of the original scheme. You always want to blend your scheme into a first order upwind (when crossing bounds)

glwagner · 2022-05-17T21:25:41Z

To support different bounds for different tracers, let's allow the advection scheme to be a NamedTuple. Then pass the advection scheme specific to each tracer into to the tracer tendency kernel. This moves us towards the design we want for #2454 and isn't any more difficult than a different solution.

simone-silvestri · 2022-05-17T21:27:51Z

I think it is already like this for the HydrostaticFreeSurfaceModel,
we can specify

c_adv = WENO5(bounds = (0, 1))
T_adv = WENO5()
model = HydrostaticFreeSurfaceModel(tracer_advection = (c = c_adv, T = T_adv), tracers = (:c, :T))

(I just found out)
I ll change the nonhydrostatic and shallow water model to have the same feature

simone-silvestri · 2022-05-17T21:42:12Z

I guess it is better to change the ShallowWaterModel in #2522 as a lot is changing already

glwagner · 2022-05-18T13:23:20Z

For NonhydrostaticModel you'll have to come up with an API / syntax for momentum versus tracers.

Definitely not a bad idea, but we don't want to abuse people with too many syntax changes so we may want to discuss what that would look like.

glwagner · 2022-05-18T13:24:06Z

I think it is already like this for the HydrostaticFreeSurfaceModel, we can specify
c_adv = WENO5(bounds = (0, 1))
T_adv = WENO5()
model = HydrostaticFreeSurfaceModel(tracer_advection = (c = c_adv, T = T_adv), tracers = (:c, :T))
(I just found out) I ll change the nonhydrostatic and shallow water model to have the same feature

Heh, I probably wrote that but then forgot! Good rediscovery. That seems like a decent API --- what do you think?

simone-silvestri · 2022-05-18T13:27:56Z

I like it! I changed the NonhydrostaticModel to be the same.

A way to facilitate users is that we leave the option to specify a global advection (so that we keep the keyword) or two separate momentum_advection and tracer_advection.
I would rather just change it to be consistent with the HydrostaticFreeSurfaceModel directly.

what do you guys think?

glwagner · 2022-05-18T13:32:48Z

I like it! I changed the NonhydrostaticModel to be the same.

A way to facilitate users is that we leave the option to specify a global advection (so that we keep the keyword) or two separate momentum_advection and tracer_advection. I would rather just change it to be consistent with the HydrostaticFreeSurfaceModel directly.

what do you guys think?

I think we should do that in a different PR because it's a major, major breaking change to the most popular model!

As for the syntax I'd prefer --- I would prefer one keyword argument advection. We need that for the changes proposed in #2454 --- because the schemes for different components can depend on one another. So it's easier to put it in one object.

I also feel it's better to just specify grid once to the advection constructor. But that would be harder to implement now.

glwagner · 2022-05-18T13:34:43Z

src/Models/NonhydrostaticModels/nonhydrostatic_model.jl

@@ -50,7 +55,8 @@ end
 """
    NonhydrostaticModel(;     grid,
                              clock = Clock{eltype(grid)}(0, 0, 1),
-                          advection = CenteredSecondOrder(),
+                 momentum_advection = CenteredSecondOrder(),
+                   tracer_advection = CenteredSecondOrder(),


I prefer one kwarg advection --- and we should do this in a different PR.

Ok, then we can change it in a different PR (for the moment then separate bounds can be used only in the HydrostaticFreeSurfaceModel)

src/Advection/positivity_preserving_tracer_advection_operators.jl

glwagner · 2022-05-20T03:02:30Z

src/Advection/positivity_preserving_tracer_advection_operators.jl

+    p̃   =  (cᵢⱼ - ω̂₁ * c₋ᴿ - ω̂ₙ * c₊ᴸ) / (1 - 2ω̂₁)
+    M   = max(p̃, c₊ᴸ, c₋ᴿ) 
+    m   = min(p̃, c₊ᴸ, c₋ᴿ) 
+    θ   = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), 1.0)


Suggested change

θ = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), 1.0)

θ = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), one(grid))

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

simone-silvestri added 2 commits May 16, 2022 17:50

first commit

b71b13a

small bug

eb5a684

simone-silvestri added 5 commits May 16, 2022 18:50

allow generic bounds

6004c61

correct gauss lobatto weights

a94d29a

tidy up a bit

8b009a0

tidy up a bit

1378102

bugfix

17cde23

added z-direction

878c015

usual Manifest/Project incompatibility...

83aef6d

simone-silvestri changed the title ~~[WIP] Maximum-principle-satisfying WENO scheme~~ Maximum-principle-satisfying WENO scheme May 17, 2022

francispoulin reviewed May 17, 2022

View reviewed changes

glwagner reviewed May 17, 2022

View reviewed changes

glwagner approved these changes May 17, 2022

View reviewed changes

francispoulin approved these changes May 17, 2022

View reviewed changes

simone-silvestri added 3 commits May 17, 2022 16:41

allow different bounds for different tracers

49595a3

allow different bounds for different tracers

9f9ce85

reverted back

2bcd03e

simone-silvestri added 2 commits May 17, 2022 17:50

momentum/tracer advection in nonhydrostatic model

d85ddbb

start bugfixing

f9f3dfc

glwagner requested changes May 18, 2022

View reviewed changes

simone-silvestri added 4 commits May 18, 2022 11:39

reverted again

2d5391c

support for flat direfctions

1c9da05

support for flat directions

63d0637

I hate mac

b020bd4

glwagner reviewed May 20, 2022

View reviewed changes

src/Advection/positivity_preserving_tracer_advection_operators.jl Outdated Show resolved Hide resolved

glwagner reviewed May 20, 2022

View reviewed changes

src/Advection/positivity_preserving_tracer_advection_operators.jl Outdated Show resolved Hide resolved

glwagner reviewed May 20, 2022

View reviewed changes

glwagner approved these changes May 20, 2022

View reviewed changes

simone-silvestri and others added 3 commits May 19, 2022 23:19

Update src/Advection/positivity_preserving_tracer_advection_operators.jl

254ed88

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

Update src/Advection/positivity_preserving_tracer_advection_operators.jl

3151ce7

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

Update positivity_preserving_tracer_advection_operators.jl

166ac86

simone-silvestri merged commit ea9a1a1 into main May 21, 2022

navidcy deleted the ss/positivity_preserving_weno branch February 8, 2023 04:01

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Maximum-principle-satisfying WENO scheme #2557

Maximum-principle-satisfying WENO scheme #2557

simone-silvestri commented May 16, 2022 •

edited

simone-silvestri commented May 16, 2022 •

edited

simone-silvestri commented May 16, 2022

francispoulin commented May 17, 2022

simone-silvestri commented May 17, 2022 •

edited

glwagner commented May 17, 2022

francispoulin May 17, 2022

simone-silvestri May 17, 2022

francispoulin May 17, 2022

glwagner May 17, 2022

glwagner May 17, 2022

simone-silvestri May 17, 2022

glwagner May 20, 2022

simone-silvestri May 20, 2022

glwagner May 20, 2022

glwagner commented May 17, 2022 •

edited

simone-silvestri commented May 17, 2022 •

edited

simone-silvestri commented May 17, 2022

glwagner commented May 18, 2022

glwagner commented May 18, 2022

simone-silvestri commented May 18, 2022 •

edited

glwagner commented May 18, 2022 •

edited

glwagner May 18, 2022 •

edited

simone-silvestri May 18, 2022

glwagner May 20, 2022

	θ = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), 1.0)
	θ = min(abs((upper_limit - cᵢⱼ)/(M - cᵢⱼ + ε₂)), abs((lower_limit - cᵢⱼ)/(m - cᵢⱼ + ε₂)), one(grid))

Maximum-principle-satisfying WENO scheme #2557

Maximum-principle-satisfying WENO scheme #2557

Conversation

simone-silvestri commented May 16, 2022 • edited

simone-silvestri commented May 16, 2022 • edited

simone-silvestri commented May 16, 2022

francispoulin commented May 17, 2022

simone-silvestri commented May 17, 2022 • edited

glwagner commented May 17, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

glwagner commented May 17, 2022 • edited

simone-silvestri commented May 17, 2022 • edited

simone-silvestri commented May 17, 2022

glwagner commented May 18, 2022

glwagner commented May 18, 2022

simone-silvestri commented May 18, 2022 • edited

glwagner commented May 18, 2022 • edited

glwagner May 18, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

simone-silvestri commented May 16, 2022 •

edited

simone-silvestri commented May 16, 2022 •

edited

simone-silvestri commented May 17, 2022 •

edited

glwagner commented May 17, 2022 •

edited

simone-silvestri commented May 17, 2022 •

edited

simone-silvestri commented May 18, 2022 •

edited

glwagner commented May 18, 2022 •

edited

glwagner May 18, 2022 •

edited