Difference between UpwindBiasedFifthOrder and WENO5 schemes

[tomchor]

This is related to a message I sent originally on Slack. Let me start by apologizing for the example being kinda long. I spent a considerable amount of time on this and I have not been able to locate the issue or come up with a shorter example. The issue also seems to be related to the part of parameter space I'm in, and related to resolution, making it hard to figure out what's happening.

Btw, I figured I'd start a discussion instead of an issue because I'm not sure it's a bug. This may be simply a numerical limitation of the schemes, or I may be missing something. Any help is appreciated regardless.

In the following example I'm simulating a 2D non-rotating domain. The whole domain is initially unstrafified but buoyancy is an active tracer. The only non-standard BC is a ValueBoundaryCondition at the west wall, which imposes a Δb (b being buoyancy) which is supposed to lead to a gravity current of sorts. The value of Δb is negative close to the bottom and zero above.

Here's the code:

using Oceananigans
using Oceananigans.Units

params = (Δb = 5e-2, # m/s²
          Δz = 40, # m
          Lx = 8e3,
          Lz = 300,
          Nx = 2^7,
          Nz = 2^3,
          )
arch = CPU()

topo = (Bounded, Flat, Bounded)
grid = RegularRectilinearGrid(topology=topo,
                              size=(params.Nx, params.Nz),
                              x=(0, params.Lx), z=(0, params.Lz),
                              halo=(3,3),
                              )
@info grid


#+++++ Buoyancy inflow
logistic(φ) = 1 / (1 + exp(-φ))
b_west_value_func(y, z, t, p) = - p.Δb * (1 - logistic((z - p.Δz)/3))
b_west_value = ValueBoundaryCondition(b_west_value_func, parameters=(Δz=params.Δz,
                                                                     Δb=params.Δb,
                                                                     ))
b_bcs = FieldBoundaryConditions(west = b_west_value,)
bcs = (b=b_bcs,)
#-----


#++++ Model and ICs
model = NonhydrostaticModel(grid = grid, timestepper = :RungeKutta3,
                            architecture = arch,
                            #advection = WENO5(),
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = BuoyancyTracer(),
                            coriolis = nothing,
                            tracers = :b,
                            closure = IsotropicDiffusivity(ν=1e-5, κ=1e-5),
                            boundary_conditions = bcs,
                           )
@info "" model
#----


#++++ Create simulation
using Oceananigans.Grids: min_Δz
cfl=0.85
wizard = TimeStepWizard(max_change=1.01, cfl=cfl)

# Print a progress message
using Oceanostics: TimedProgressMessenger
simulation = Simulation(model, Δt=min(1, 0.5*min_Δz(grid)/1e-3), stop_time=2days,)
simulation.callbacks[:progress] = Callback(TimedProgressMessenger(LES=false,), IterationInterval(10))
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(2))
#----


#++++ Diagnostics
u, v, w = model.velocities
b = model.tracers.b
outputs = (u=u, v=v, w=w, b=b,)
simulation.output_writers[:vid_writer] =
    NetCDFOutputWriter(model, outputs,
                       filepath = "data/vid.onefile.nc",
                       schedule = TimeInterval(20minutes),
                       mode = "c",
                       )
#----


@info "Starting simulation"
run!(simulation, pickup=false)

Now, using the Upwind scheme, this is what the solution looks like for the buoyancy b:

b_upw5.mp4

Now, because I was using this for many tests it's pretty under-resolved, which is probably what's making the circulation kinda weird, but at least I see a circulation develop, which is what I expected. (Although I don't know why there's apparently both "heating" above and "cooling" below, since Δb is non-positive...)

If I run the exact same example with a WENO5 advection scheme, this is what I get:

b_weno5.mp4

In order words: I get no circulation at all. As if all BCs were no flux.

This happened originally with a code of mine that's more complicated and which is better resolved, but the gist is the same. Am I missing something?

CC @glwagner @navidcy

[navidcy]

Can you post the plotting code. I can try to reduce this to something more minimal. That'll help us.

E.g., do you even need non-flat z?

[tomchor]

I actually found that, for this set of parameters, a non-tilted domain works. Also no need for stratification I guess. So I updated the example to something a bit more minimal.

So in the bottom WENO5 case, the model.u velocity remains zero at all time? Is that the "issue"?

Yes, the issue is that with WENO5 the flow stays quiescent despite the forcing on the west wall. (All the velocities are approximately zero.)

[tomchor]

Also, just to make sure, this is with a tagged Oceananigans version? Which one?

Another good point, here's my ]st:

(dcurrent4) pkg> st
      Status `/glade/scratch/tomasc/dcurrent4/Project.toml`
  [c7e460c6] ArgParse v1.1.4
  [052768ef] CUDA v3.3.6
  [85f8d34a] NCDatasets v0.11.7
  [9e8cae18] Oceananigans v0.65.0
  [d0ccf422] Oceanostics v0.6.1 `https://github.com/tomchor/Oceanostics.jl.git#main`
  [d96e819e] Parameters v0.12.3
  [91a5bcdd] Plots v1.23.5
  [d330b81b] PyPlot v2.10.0

[navidcy]

Yes, the issue is that with WENO5 the flow stays quiescent despite the forcing on the west wall. (All the velocities are approximately zero.)

approximately or exactly zero?

[tomchor]

Approximately. Max velocities for the Upwind scheme are about 7e-1 m/s, while for WENO5 they're about 3e-2 m/s.

[glwagner]

It's possible that the WENO5() result more accurately reflects the true solution.

[glwagner]

When using ValueBoundaryCondition, any boundary fluxes that occur are tied to the specified turbulence closure. In this example, we'll see fluxes into the domain because of the inhomogeneous ValueBoundaryCondition combined with closure = IsotropicDiffusivity(ν=1e-5, κ=1e-5).

My best guess is that the behavior observed when using advection=UpwindBiasedFifthOrder() is associated with numerical errors. WENO5() is more accurate and tends to produce future spurious oscillations; therefore the solution when using WENO5() may be more reliable. One indication this could be the case is that it appears that the net flux is similar in both simulations (eg in the UpwindBiasedFifthOrder() case, we see equal amounts of positive and negative buoyancy, so the total flux is small). It's also suspicious that the diffusivity is quite small at 1e-5. How much buoyancy flux do we expect in this case given the other parameters in the problem? I think changing the colorscale for the WENO5() case could confirm that there is indeed buoyancy flux (because of the non-zero diffusivity and ValueBoundaryCondition), but that the buoyancy perturbations associated with this flux are too small to induce any noticeable circulation (eg, the Peclet number is small).

[navidcy]

To remedy all these I'd feel better if there was some analytic solution for the particular boundary condition we can then reproduce with the two advection schemes. A one-dimensional one x with no-flux BC on right and prescribed value on the left. And with a constant u... Just brainstorming here.

Otherwise we'd be in constant limbo... My 2 cents.

[tomchor]

I agree that it looks like numerical errors, and I actually considered that, but aren't all boundary fluxes implemented with the same scheme? I remember Ali (and I think you too) mentioning that value boundary conditions were kinda tricky, so we limited to implementing them with a simpler operator always? (Or maybe I'm confusing with advection near the wall...)

[glwagner]

I agree that it looks like numerical errors, and I actually considered that, but aren't all boundary fluxes implemented with the same scheme? I remember Ali (and I think you too) mentioning that value boundary conditions were kinda tricky, so we limited to implementing them with a simpler operator always? (Or maybe I'm confusing with advection near the wall...)

When using ValueBoundaryCondition, any boundary fluxes that occur are tied to the specified turbulence closure; the turbulence closure computes the flux.

[glwagner]

As a side note, if indeed UpwindBiasedFifthOrder() produces such radically lower quality results as WENO5() (qualitatively different!), this is perhaps more fuel to make WENO5() the default advection scheme! If only it weren't so slow on the CPU...

[navidcy]

And imagine that MOM5 uses 2nd order advection (I think).

[tomchor]

this is perhaps more fuel to make WENO5() the default advection scheme! If only it weren't so slow on the CPU...

Is it worth implementing a WENO3 scheme? It might be a good compromise between cost and accuracy. I'm assuming it wouldn't be hard to implement, since we already have WENO5.

[navidcy]

OK, I can reproduce this.

But my question is why one is more correct than the other? I have no reason to believe that either is more correct than the other. They are indeed different, but should they be same?

That's why I'd feel that if we don't have a sense of how the solution should be then why go after chasing a bug in either of the two? @tomchor you seem certain that a bug, if any, should be in WENO5. Perhaps you have more intel. But just from these two exps I'm not sure why it's not the UpwindBiasedFifthOrder scheme that is at fault.

Both the animations show buoyancy with limits (-1e-3, 1e-3).

using Oceananigans
using Oceananigans.Units

advection_scemes = (WENO5(), UpwindBiasedFifthOrder())

for advection_sceme in advection_scemes

if advection_sceme == WENO5()
    filename = "test_tomas_WENO5"
elseif advection_sceme == UpwindBiasedFifthOrder()
    filename = "test_tomas_Upwind5"
end

params = (Δb = 5e-2, # m/s²
          Δz = 40, # m
          Lx = 8e3,
          Lz = 300,
          Nx = 2^7,
          Nz = 2^3,
          )

arch = CPU()

topo = (Bounded, Flat, Bounded)

grid = RegularRectilinearGrid(topology=topo,
                              size=(params.Nx, params.Nz),
                              x=(0, params.Lx), z=(0, params.Lz),
                              halo=(3,3))

logistic(φ) = 1 / (1 + exp(-φ))

b_west_value_func(y, z, t, p) = - p.Δb * (1 - logistic((z - p.Δz)/3))
b_west_value = ValueBoundaryCondition(b_west_value_func, parameters=(Δz=params.Δz,
                                                                     Δb=params.Δb,
                                                                     ))
b_bcs = FieldBoundaryConditions(west = b_west_value,)
bcs = (b=b_bcs,)

model = NonhydrostaticModel(grid = grid, timestepper = :RungeKutta3,
                            architecture = arch,
                            advection = advection_sceme,
                            buoyancy = BuoyancyTracer(),
                            coriolis = nothing,
                            tracers = :b,
                            closure = IsotropicDiffusivity(ν=1e-1, κ=1e-1),
                            boundary_conditions = bcs)

simulation = Simulation(model, Δt=5, stop_time=1days)

u, v, w = model.velocities
b = model.tracers.b

outputs = (u=u, v=v, w=w, b=b,)

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                      schedule = TimeInterval(20minutes),
                                                      prefix = filename,
                                                      field_slicer = nothing,
                                                      verbose = false,
                                                      force = true)

@info "Starting simulation"
run!(simulation, pickup=false)

using GLMakie

filepath = filename * ".jld2"

bt = FieldTimeSeries(filepath, "b")

x, y, z = nodes((Center, Center, Center), grid)

times = bt.times
Nt = length(times)

bn(n) = interior(bt[n])[:, 1, :]

n = Node(1)

b = @lift bn($n)

fig = Figure()
ax = Axis(fig[1, 1], title="buoyancy")
hm = heatmap!(ax, x, z, b, colorrange=(-1e-3, 1e-3), colormap=:balance)

display(fig)

record(fig, filename * ".mp4", 1:Nt, framerate=8) do i
    @info "Plotting frame $i of $Nt"
    n[] = i
end

end

(tomas_debug) pkg> st
      Status `~/Research/tomas_debug/Project.toml`
  [e9467ef8] GLMakie v0.4.7
  [9e8cae18] Oceananigans v0.65.0

WENO5

test_tomas_WENO5.mp4

UpWind5

test_tomas_Upwind5.mp4

[glwagner]

(Almost two orders of magnitude difference in the velocities.)

We don't have to look far to come up with a problem that blows up with low-order advection, but runs stably with high-order advection. Such a difference is far greater than a few orders of magnitude!

An issue with all finite volume schemes is that numerical errors can corrupt the solution even when the time-integration doesn't produce NaNs. This is in contrast to a spectral method that might simply blow up when a physical/numerical simulation is setup poorly (for example by under-resolving dissipative scales, or using too-large a time step). I think we should be pleased that WENO5 appears to produce a sensible result even at such a coarse resolution.

Could the 2nd order advection near the wall be influencing this? @glwagner can you explain why we do that? This definitely introduces errors in solutions.

All advection schemes limit to 2nd order near boundaries (and all schemes introduce errors, not just 2nd order ones). Since the limiting is the same for both advection schemes here, I don't see why this specifically causes the UpwindBiasedFifthOrder case to become corrupted.

Here's a brief explanation of the finite volume numerical method being used for high-order reconstructions in a Bounded direction:

2nd order reconstruction is the only reconstruction we have currently implemented that's valid in boundary-adjacent cells.

It is possible to implement a different algorithm that replaces an nth order "interior reconstruction" with an nth order "inward looking" reconstruction near the boundary. This algorithm would maintain the order of accuracy of the advective reconstruction near the boundary and would be a great contribution to Oceananigans. Hopefully it would also improve the fidelity of solutions near the boundary. The reason I'm not sure it would help is that 1) the reconstruction would be "downwind" (when normal velocities expel fluid from a boundary) as often as it would be "upwind" (when normal velocities impinge on a boundary), which might destabilize the solution. 2) advective fluxes normal to boundaries are often smaller than tangential advective fluxes (which are not affected by this discussion); therefore overall error may be controlled by tangential reconstruction near boundaries.

There could be other algorithms worth entertaining too. For example another possibility is to fill halo regions as many times as needed for different reconstructions (there's some arithmetic required to figure out how to fill halo regions correctly for higher-than-second-order stencils), and then to evaluate terms with different stencils in the tendency in separate kernels. But this be complicated to implement --- and I'm not sure this algorithm actually yields a higher-than-second-order numerical scheme.

[tomchor]

An issue with all finite volume schemes is that numerical errors can corrupt the solution even when the time-integration doesn't produce NaNs. This is in contrast to a spectral method that might simply blow up when a physical/numerical simulation is setup poorly (for example by under-resolving dissipative scales, or using too-large a time step). I think we should be pleased that WENO5 appears to produce a sensible result even at such a coarse resolution.

I guess that's a good point. Most of my intuition comes from spectral methods, but things are different for finite volumes. One way to investigate this further (for anyone who's got the resources) is to refine both solutions. If the issue is, indeed, numerics, they should both converge, no?

[glwagner]

If the issue is, indeed, numerics, they should both converge, no?

Yes I think so. But here the Peclet / Reynolds number is around 10^5 (U=1e-2, L=300, kappa=1e-5). So we might be a long ways from resolving the dissipation scale here and achieving numerical convergence?

[navidcy]

@glwagner, note that I used different κ, v values for the movies I made. But perhaps that's why they look more similar to the original that @tomchor posted?

[glwagner]

Oh true. Yeah, Re = Pe = 1000 is doable.