bug with background fields?

[johnryantaylor]

Hi all, We have come across what seems to be a bug in background_fields. This arose in another context, but the code below is the simplest way that I have come up with so far that demonstrates the problem. Here there is a stable linear buoyancy profile (imposed as a background gradient) with random noise for the initial velocity. This generates extra turbulence at the top and bottom. I wonder if the background gradient in the halo region isn't playing nicely with the initial conditions. Has anyone else seen something like this or have any ideas for where to look for the problem? I am using Julia 1.9.2 and Oceananigans 0.85.0 Thanks! John

using Oceananigans
using JLD2
using Plots
using Printf

Lx = 1
Lz = 1
Nx = 64
Nz = 64

grid = RectilinearGrid(size=(Nx, Nz), x=(0, Lx), z=(0, Lz), topology=(Periodic, Flat, Periodic))

# Create a stable background buoyancy gradient
N = 1       # buoyancy frequency [s⁻¹]
B_func(x, y, z, t, N) = N^2 * z
B = BackgroundField(B_func, parameters=N)

model = NonhydrostaticModel(; grid,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            closure = ScalarDiffusivity(ν=1e-6, κ=1e-6),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            background_fields = (; b=B)) # `background_fields` is a `NamedTuple`

A=0.01
u_ic(x,y,z) = A * randn()
w_ic(x,y,z) = A * randn()

set!(model, u=u_ic, w=w_ic)

simulation = Simulation(model, Δt = 0.01, stop_iteration=200)

wizard = TimeStepWizard(cfl=0.3, diffusive_cfl=0.3, max_change=1.2, max_Δt=1)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))

progress(sim) = @printf("i: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, advective CFL: %.2e, diffusive CFL: %.2e\n",
                        iteration(sim), time(sim), prettytime(sim.run_wall_time),
                        sim.Δt, AdvectiveCFL(sim.Δt)(sim.model), DiffusiveCFL(sim.Δt)(sim.model))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

filename = "background_field_test"
simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          filename = filename * ".jld2",
                                                          schedule = IterationInterval(10),
                                                          overwrite_existing = true)

run!(simulation)

xu, yu, zu = nodes(model.velocities.u) 
xw, yw, zw = nodes(model.velocities.w)

file = jldopen(filename * ".jld2")

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i,iter) in enumerate(iterations)
    @info "Drawing frame $i from iteration $iter..."

    u = file["timeseries/u/$iter"][:,1,:];
    w = file["timeseries/w/$iter"][:,1,:];
    t = file["timeseries/t/$iter"];

    ## Plot the heatmaps
    u_plot = Plots.heatmap(xu, zu, u', color=:thermal, xlabel="x", ylabel="z", aspect_ratio = :equal, xlims = (0, Lx), ylims = (0, Lz));
    w_plot = Plots.heatmap(xw, zw, w', color=:thermal, xlabel="x", ylabel="z", aspect_ratio = :equal, xlims = (0, Lz), ylims = (0, Lz));
    u_title = @sprintf("u, t=%s",round(t));
    w_title = @sprintf("w, t=%s",round(t));

    Plots.plot(u_plot,w_plot,layout=(1,2),size=(1200,600),title =[u_title w_title]);

    iter == iterations[end] && close(file) # close the file if we reach the end

    end

mp4(anim, filename * ".mp4", fps = 10);

[glwagner]

Hmm, I think there could be a bug with vertically-periodic simulations, since here topology = (Periodic, Flat, Periodic). @tomchor mentioned finding something similar I think? I believe we have to fix something with the way the hydrostatic pressure component is treated (or eliminate it completely).

I'm not quite sure how issues with being vertically-periodic interact with the background field here but there might be something...

[glwagner]

It could be worth trying to run on this branch https://github.com/CliMA/Oceananigans.jl/pull/3080

[tomchor]

This is the issue that described what I was seeing: https://github.com/CliMA/Oceananigans.jl/issues/3290

@johnryantaylor does that looks like what you're seeing?

If so, the two workarounds that I have come up with are

• to run on branch #3080 like @glwagner suggested (downside here is that this branch is kinda out-of-date and apparently stale)
• Make the x (or y) directions periodic (seems like xalready is in your case) and tilt the domain accordingly. (The second video in https://github.com/CliMA/Oceananigans.jl/issues/3290 illustrates that)

If that's not your issue, do you mind posting the animatino you're generating?

[glwagner]

Given that there seems to be increasing interest in vertically-periodic simulations, we could revive #3080 (despite the caveats mentioned there, which mainly unknowns associated with the performance of potential future nonhydrostatic solvers for complex domains), since it's always possible to reverse course in the future and restore the separation (perhaps when the separation is restored, it can be done in a way that's compatible with vertically periodic domains).

[glwagner]

This is the issue that described what I was seeing: #3290

@johnryantaylor does that looks like what you're seeing?

If so, the two workarounds that I have come up with are

• to run on branch Don't separate the pressure into hydrostatic and nonhydrostatic in NonhydrostaticModel #3080 like @glwagner suggested (downside here is that this branch is kinda out-of-date and apparently stale)
• Make the x (or y) directions periodic (seems like xalready is in your case) and tilt the domain accordingly. (The second video in Having a domain that's periodic in the "gravity" direction doesn't work for gravity_unit_vector = NegativeZDirection() #3290 illustrates that)

If that's not your issue, do you mind posting the animatino you're generating?

And just to confirm @tomchor those issues are not associated with BackgroundField, right?

[glwagner]

I can confirm no issues with vertically Bounded.

Here's slightly modified code

using Oceananigans
using GLMakie
using Printf

grid = RectilinearGrid(size = (64, 64),
                       x = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Flat, Periodic))

@inline B_func(x, y, z, t) = z
B = BackgroundField(B_func)

model = NonhydrostaticModel(; grid,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            closure = ScalarDiffusivity(ν=1e-6, κ=1e-6),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            background_fields = (; b=B)) # `background_fields` is a `NamedTuple`

A = 0.01
u_ic(x, y, z) = A * randn()
w_ic(x, y, z) = A * randn()

set!(model, u=u_ic, w=w_ic)

simulation = Simulation(model, Δt=0.01, stop_iteration=200)

wizard = TimeStepWizard(cfl=0.3, diffusive_cfl=0.3, max_change=1.2, max_Δt=1)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))

progress(sim) = @printf("i: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, advective CFL: %.2e, diffusive CFL: %.2e\n",
                        iteration(sim), time(sim), prettytime(sim.run_wall_time),
                        sim.Δt,
                        AdvectiveCFL(sim.Δt)(sim.model),
                        DiffusiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

filename = "background_field_test.jld2"
simulation.output_writers[:velocities] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          filename = filename,
                                                          schedule = IterationInterval(10),
                                                          overwrite_existing = true)

run!(simulation)

wt = FieldTimeSeries(filename, "w")
bt = FieldTimeSeries(filename, "b")
t = wt.times
Nt = length(t)
n = Observable(1)

wn = @lift interior(wt[$n], :, 1, :)
bn = @lift interior(bt[$n], :, 1, :)

fig = Figure(resolution=(1500, 700))
axw = Axis(fig[1, 1])
axb = Axis(fig[1, 2])

heatmap!(axw, wn)
heatmap!(axb, bn)

topostr = string(Oceananigans.Grids.topology(grid, 3))

record(fig, "vertically_$topostr.mp4", 1:Nt, framerate=12) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end

vertically_Bounded.mp4

https://github.com/CliMA/Oceananigans.jl/assets/15271942/8b63ae62-1a44-47c2-931b-9db154adddc2

vertically_Periodic.mp4

https://github.com/CliMA/Oceananigans.jl/assets/15271942/8af44a0f-451e-4f58-a6f0-e204a4e65bd3

The easiest fix is to eliminate the pressure separation.

Another solution is to fix the hydrostatic pressure algorithm.

We should also note that the vertical tridiagonal solve is not correct for vertically-periodic domains. But this is easily solvable (and the above 2 are as well).

[tomchor]

This is the issue that described what I was seeing: #3290 @johnryantaylor does that looks like what you're seeing? If so, the two workarounds that I have come up with are

• to run on branch Don't separate the pressure into hydrostatic and nonhydrostatic in NonhydrostaticModel #3080 like @glwagner suggested (downside here is that this branch is kinda out-of-date and apparently stale)
• Make the x (or y) directions periodic (seems like xalready is in your case) and tilt the domain accordingly. (The second video in Having a domain that's periodic in the "gravity" direction doesn't work for gravity_unit_vector = NegativeZDirection() #3290 illustrates that)

If that's not your issue, do you mind posting the animatino you're generating?

And just to confirm @tomchor those issues are not associated with BackgroundField, right?

If it's the same issue as the one I pointed out in #3290 (which is not 100% clear atm), then you're right, they're not associated with BackgroundField. They have to do with trying to enforce a periodic pressure when the hydrostatic pressure is aperiodic on z due to the fact that it comes from a b integral in z.

[glwagner]

For ocean applications this would probably be the most important application (ie interior ocean turbulence in the presence of a buoyancy gradient, which can be modeled as vertically-periodic if we simulate perturbations about the background gradient).

[tomchor]

Thanks for running these tests, @glwagner. They really do look like the issue in https://github.com/CliMA/Oceananigans.jl/issues/3290

The easiest fix is to eliminate the pressure separation.

Since we're waiting for the IBM-aware pressure solve, we could also just add an option to eliminate the pressure separation that would be false by default. That should be pretty easy. And then in the future when we're confident about the new algorithm we can eliminate the pressure separation completely (along with user interface an code simplifications that are possible with no pressure separation).

Another solution is to fix the hydrostatic pressure algorithm.

Does it need fixing in this case though? The way I see it this is just a consequence of how the hydrostatic pressure is defined: a vertical integral of b, which doesn't play well with the assumption of a vertically-periodic domain.

We should also note that the vertical tridiagonal solve is not correct for vertically-periodic domains. But this is easily solvable (and the above 2 are as well).

Cool! We should probably do that as well :)

[glwagner]

we could also just add an option to eliminate the pressure separation that would be false by default.

I think we already have a keyword argument for pressure_solver, so we can generalize that to indicate separation vs not. I would also change "pressures" property to just pressure (which can still be NamedTuple if we are using a separation).

[glwagner]

Does it need fixing in this case though? The way I see it this is just a consequence of how the hydrostatic pressure is defined: a vertical integral of b, which doesn't play well with the assumption of a vertically-periodic domain.

I don't think there is anything essentially incompatible about using an integral per se. The problem is that the current implementation is not compatible with vertical periodicity. Just like the vertical tridiagonal solver, the boundary conditionn / topology is currently implicitly incorporated into the numerical method used for the integral.

[johnryantaylor]

Hi both, Thanks for the quick replies on this. I agree that it seems like it is likely coming from the background gradient and lack of periodicity when applied to the hydrostatic presssure. As far as I can tell there aren't any problems for passive scalars which supports that interpretation. As an aside, I think that there are also boundary artifacts in the internal_wave.jl example (although you need to run the simulation longer and you can't see them with the default contouring of the plots).

A solution that should work is to implicitly cancel the hydrostatic pressure gradient associated with the background field with the buoyancy term. This is what we do in Diablo (although we don't decompose the pressure into hydrostatic and nonhydrostatic components). Implementing this is simple since you just don't include the background gradient when calculating the hydrostatic pressure. The only trouble that I see is figuring out when to do this. In other contexts you want to keep this term. E.g. in the geostrophic adjustment problem, the hydrostatic pressure gradient is needed to drive the flow. Maybe it would be safe to exclude the background buoyancy field in the calculation of the hydrostatic pressure anytime when the topology is periodic in the vertical direction?

[johnryantaylor]

Actually, I think I misunderstood the problem... Based on update_hydrostatic_pressure.jl, it doesn't look like the background buoyancy field is used in the calculation of the hydrostatic pressure. Indeed a quick test with a background horizontal buoyancy gradient and no initial velocity doesn't drive a flow. Now I think I understand what @tomchor and @glwagner were saying: the problem arises due to the vertical periodicity and not the background buoyancy field.

[glwagner]

Just for completeness I'll respond to the first comment, since I think the way we include background fields deserves more discussion (I'm not sure it's what we always want to do)

A solution that should work is to implicitly cancel the hydrostatic pressure gradient associated with the background field with the buoyancy term. This is what we do in Diablo (although we don't decompose the pressure into hydrostatic and nonhydrostatic components).

We prescribe the background fields to enter through advection terms (essentially assuming that advection is the "key nonlinearity", in the sense outlined in the docs equations). This might not be what one wants for some terms though, for example, one might want to include background shear in their LES viscosity calculation. (Possibly it would not be difficult to include background shear in a viscosity calculation, because now we have the total_velocities / SumOfArrays abstraction).