Initially quasi-quiescent stably stratified simulations develop instabilies and eventually blow up

[tomchor]

I'm seeing a pattern when running stratified simulations where I get (numerical) instabilities that keep amplifying until they eventually blow up. This can be illustrated with the following example:

using Oceananigans
using Oceananigans.Units

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded), size = (64, 32), x = (0, 400), z = (0, 100))

model = NonhydrostaticModel(; grid, timestepper = :RungeKutta3,
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b,
                            )

noise(x, y, z) = 1e-6 * rand()
@inline constant_stratification(x, y, z) = 1e-5 * z
set!(model, u = noise, b = constant_stratification)

using Statistics: mean
model.velocities.u.data.parent .-= mean(model.velocities.u)

simulation = Simulation(model, Δt = 100, stop_time = 1day)
wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))

B = Field(model.tracers.b + model.background_fields.tracers.b)
simulation.output_writers[:fields] = NetCDFOutputWriter(model, (; model.velocities..., B);
                                                        filename = "bbl.nc",
                                                        schedule = TimeInterval(10minutes),
                                                        overwrite_existing = true)
run!(simulation)

The example above sets up a 2D (xz) domain in which b is initially (stably) linearly stratified . All the boundary conditions are default, and the only other IC modification is to add a zero-mean small-scale noise (1e-6 m/s) to u. I'd expect this to remain quasi-quiescent. Instead, this is what I get:

bbl.nc.mp4

The code to plot the animation can be found here.

The same thing happens when b=0 initially and the linear stratification is instead included as a BackgroundVelocity. For example, if I remove the drag boundary condition from the tilted bottom boundary layer example and add the same zero-mean noise, the same issue happens:

tilted_bbl.nc.mp4

Here's the code for the example above.

I understand that hydrostatic balance in a mathematical sense doesn't exactly translate directly to a discrete system (since the forces might not exactly cancel due to the discretization, interpolation, etc.), but I also don't expect deviations from it that are this large and keep amplifying. In fact, what led me here was some pretty strange behavior that I'm experiencing in one of my simulations where internal waves seem to appear out of nowhere, so I believe this is relevant.

Am I missing something here? Or is this a bug?

An important note is that we used to have a test that spun up a simulation much like the one above, ran it for a few hours, and made sure that there wasn't any movement. That test used to pass but was remove because it was too computationally-expensive. So I think this issue might have been introduced relatively recently.

EDIT:

Note that this happens even when the added noise is (on average) well below the machine epsilon (for example for noise(x, y, z) = 1e-18 * rand()).

[navidcy]

But it's not initially quiescent, there is some flow.

[tomchor]

But it's not initially quiescent, there is some flow.

That's a fair point. I updated the issue to be more precise. I also added a note at the bottom to report that this happens even when the average noise is below machine-precision.

In any case, I don't this behavior of amplifying noise is to be expected regardless of the noise amplitude, no? I should also say that this happens with other time-steppers and advection schemes (including high order WENO schemes).

[glwagner]

Is energy increasing? If so that might be a time stepping error and you need to reduce your time step.

[tomchor]

Is energy increasing? If so that might be a time stepping error and you need to reduce your time step.

Yes, and the issue does seem to decrease with a time step but not nearly fast enough (at least for this specific setup). Here are some results with the MWE above using 1e-6 m/s initial noise amplitude and running the simulation for 10 days.

In all cases the initial KE is ~5e-14 m2/s2.

Target CFL	KE final / KE initial
5e-1	1e7
5e-2	1e5
5e-3	1e3
5e-4	16

So even running with a 5e-4 cfl there's still an order of magnitude increase in the KE after 10 days. To me this is totally unexpected for the setup in question.

[glwagner]

What is the source of energy? (Of course it's unexpected for a simulation to blow up --- unless you like simulating NaNs??)

[tomchor]

What is the source of energy? (Of course it's unexpected for a simulation to blow up --- unless you like simulating NaNs??)

In this case, unless I'm missing something, there shouldn't be a source of energy. All the BCs are no-flux and there are no forcing terms. So the only thing I can think of is a numerical issue.

[glwagner]

Have you tried fixing the time-step and adjusting that (rather than setting target CFL)?

[glwagner]

Note the buoyancy oscillation time-scale here is 316 s

julia> N² = 1e-5
1.0e-5

julia> 1 / sqrt(N²)
316.2277660168379

so I think that the initial dt=100 is a little close to that. You may want to try dt=10. Also, the wave time-scale is not captured by CFL / time step wizard --- you have to restrict the time step also in the wizard with something like

wizard = TimeStepWizard(max_change=1.1, cfl=0.7, max_Δt = 0.1 / sqrt(N²))

[tomchor]

That's a good point. Upon further inspection I can see that indeed the time step wizard increases Δt past that 100 s at some point.

I ran the same set-up for a fixed Δt=10 (also for 10 days) and the energy increased by a factor of 5. Definitely better! But it still seems too large of an increase for a quasi-quiescent run for such a short period of time for me. I'd expect something in the single-digit percentages...

[tomchor]

For future reference, this is what the KE evolution looks like for the simulation with Δt=10 s:

And, for comparison, this is what is looks like for the same simulation but with buoyancy=nothing:

The latter figure looks way more reasonable to me (with only about 2.5% of energy lost due to numerics).

[glwagner]

How could energy increase at all? This is stably stratified right?

[tomchor]

How could energy increase at all? This is stably stratified right?

Yup, initial distribution of buoyancy is pretty much just 1e-5 * z.

Like I said, I haven't been able to figure out why energy is increasing. Although it's hard for me to imagine anything other than a numerical instability.

[navidcy]

How could energy increase at all? This is stably stratified right?

I guess that's the whole point of the issue?

[navidcy]

For future reference, this is what the KE evolution looks like for the simulation with Δt=10 s:

And, for comparison, this is what is looks like for the same simulation but with buoyancy=nothing:

The latter figure looks way more reasonable to me (with only about 2.5% of energy lost due to numerics).

At the top figure, Initially it falls! Then it starts going up. What's happening there? That might be illuminating.

[navidcy]

You start with a random flow which is not incompressible and I presume it becomes incompressible after the first time step.

If you start with non-zero flow, why do you expect that flow to die out?

[glwagner]

Numerical dissipation by UpwindBiasedFifthOrder() should dissipate energy, similarly as in the two dimensional turbulence example.

I would also expect to see a lot of small scale waves because this is stably stratified. I guess those would also dissipate?

[glwagner]

Here's a script with a visualization

using Oceananigans
using Oceananigans.Units
using GLMakie

grid = RectilinearGrid(topology=(Periodic, Flat, Bounded), size=(64, 32), x=(0, 200), z=(0, 100))

model = NonhydrostaticModel(; grid,
                            timestepper = :RungeKutta3,
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b)

noise(x, y, z) = 1e-6 * (rand() - 1/2)
constant_stratification(x, y, z) = 1e-5 * z
set!(model, u=noise, b=constant_stratification)

simulation = Simulation(model, Δt=10.0, stop_iteration=1000)
progress(sim) = @info string("Iter: ", iteration(sim), ", time: ", prettytime(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

u, v, w = model.velocities
e = (u^2 + w^2) / 2
E = Field(Integral(e, dims=:))
outputs = merge(model.velocities, model.tracers, (; e, E))

filename = "test.jld2"
simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs;
                                                      filename,
                                                      schedule = IterationInterval(10),
                                                      overwrite_existing = true)

run!(simulation)

wt = FieldTimeSeries(filename, "w")
et = FieldTimeSeries(filename, "e")
bt = FieldTimeSeries(filename, "b")
Et = FieldTimeSeries(filename, "E")
times = wt.times
Nt = length(times)

fig = Figure(resolution=(1800, 900))

axw = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)", title="Vertical velocity")
axe = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)", title="Kinetic energy")
axb = Axis(fig[2, 3], xlabel="x (m)", ylabel="z (m)", title="Buoyancy")
axE = Axis(fig[3, 1:3], xlabel="Time (s)", ylabel="Volume averaged kinetic energy")

slider = Slider(fig[1, 1:2], range=1:Nt, startvalue=1)
n = slider.value

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

x, y, z = nodes(et)

wlim = maximum(abs, wt)
elim = maximum(abs, et)

heatmap!(axw, x, z, wn, colorrange=(-wlim, wlim), colormap=:balance)
heatmap!(axe, x, z, en, colorrange=(0, elim), colormap=:solar)
heatmap!(axb, x, z, bn, colorrange=(0, 1e-3), colormap=:thermal)

t = @lift times[$n]
lines!(axE, times, Et[:])
vlines!(axE, t)

display(fig)

record(fig, "test.mp4", 1:Nt) do nn
    @info "Recording frame $nn of $Nt"
    n[] = nn
end

github won't let me attach the movie or image but basically at first energy decreases, then there is a lot of grid scale noise with wave dynamics. So nothing is resolved, and the energy oscillates around. Perhaps later it increases. This might need a background diffusivity?

[navidcy]

If model is prescribed with just buoyancy = BuoyancyTracer() then the diffusivity is assumed zero?

[glwagner]

The default is closure=nothing, regardless of buoyancy.

[glwagner]

But there's a big difference in the dynamics of course; with buoyancy=nothing then we just have 2D turbulence whereas with stratification (very strong here relative to the initial velocity perturbations) the energy does not dissipate as strongly.

[navidcy]

test.mp4

[glwagner]

I guess when energy decreases to the point that advection is negligible compared to buoyancy/wave dynamics, the numerical diffusivity has no effect. At that point we just have wave dynamics acting on grid scale noise. I don't know what we expect. It is surprising that things blow up eventually. Adding closure = ScalarDiffusivity(ν=1e-4, κ=1e-4) causes energy to decrease monotonically.

[glwagner]

Here's what happens if you put in background diffusivity

test.mp4

[glwagner]

The buoyancy is actually a buoyancy perturbation, showing the growth of diffusive boundary layers + waves. But the viscosity and diffusivity help dissipate the grid scale waves so that larger scale waves dominate later on.

[glwagner]

Should we convert this to discussion? I don't think there's a bug in the code, but it does show that you may need background diffusivity for weakly nonlinear simulations dominated by waves.

[glwagner]

With no background diffusivity I also find the energy increases slowly. This is with WENO advection scheme.

test.mp4

Perhaps one could compute the discrete conservation laws for the linear wave equations here given our staggered discretization of buoyancy / hydrostatic pressure + nonhydrostatic pressure... I'm not sure what we would find. It's not that simple, because nonhydrostatic pressure is treated implicitly / with a fractional step. Our hydrostatic pressure scheme is not the only possible scheme one might use. Perhaps other schemes have different numerical properties.

[glwagner]

Another solution is to increase the amount of energy in the initial condition, which leads to strong numerical dissipation at the beginning and yields a smoother solution at later times...

[tomchor]

At the top figure, Initially it falls! Then it starts going up. What's happening there? That might be illuminating.

Numerical dissipation by UpwindBiasedFifthOrder() should dissipate energy, similarly as in the two dimensional turbulence example.

I would also expect to see a lot of small scale waves because this is stably stratified. I guess those would also dissipate?

I think the initial big drop has to do with energy being quickly exchange between KE and APE in this case. So I think that big drop is expected given the initial small-scale noise we include in the IC.

[tomchor]

Should we convert this to discussion? I don't think there's a bug in the code, but it does show that you may need background diffusivity for weakly nonlinear simulations dominated by waves.

I'm not sure I can rule out a bug in the code. Can you confirm that you ran a case where you get a final KE that's (significantly) larger than the initial KE? If so, then I think it's hard to imagine a situation where that would happen without some sort of bug in the code.

If the final KE is approximately equal or smaller than the initial KE, then I agree this might just be a dynamical behavior that's kinda unexpected.

[glwagner]

What's your reasoning that there is a bug in the code? Can you prove that KE won't increase with a finite time-step given the second-order discretization + fractional step for pressure?

[glwagner]

In this case if we see convergence with decreasing time-step (spatial resolution isn't relevant for these dynamics, which are at the grid scale and therefore not physical) then I suppose that would be an indication the issue is due to a finite time step.

It's important to recognize that the dynamics of the discrete linear equations are different than the continuous. When we have a smooth solution, such that our spatial discretization should approximate some exact smooth solution, we can test that refining the grid and time step leads to convergence to an exact solution. Moreover, we can estimate the time-scale of the dynamics using the time-scales of the smooth dynamics as a guide. An example is a resolved buoyancy oscillation: it has a timescale of roughly 1/N.

This example is dominated by small amplitude (eg linear) noise at the grid scale. Therefore my initial time scale estimate of 1/N may not hold. Instead, we'd have to look at the discrete eigenvalues of the system at the very higheset wavenumbers (ie the Nyquist number 2pi / dx). We could then calculate the time-step that would be required to resolve these (completely unphysical) dynamics.

I'm not sure what a lower bound on such spurious discrete dynamics might be. It could be far smaller than any physical time scale, ie as small as 1e-16, or smaller? Perhaps the evolution of grid scale noise also has to do with spatial resolution, so that could be another knob to vary.

So if we want to investigate this further, we should conduct a systematic study of the dynamics of this grid-scale noise system affected by buoyancy, decreasing the time-step to zero.

Or we can convince ourselves that non-noisy dynamics are accurate --- eg by analyzing a system like the one we use for our "internal wave" dynamics test:

https://github.com/CliMA/Oceananigans.jl/blob/main/test/test_internal_wave_dynamics.jl

that test verifies that a wave packet in our code propagates at the correct group speed, for example. But one could dive much deeper into that example and test a wide variety of wave numbers.

It depends whether you are interested in the grid-scale noise system, or whether you are interested in verifying that smooth resolved dynamics are correct.

[navidcy]

I agree. I would be more alarmed if an initial condition that includes small scales but resolved (that is with 5-6 grid points per shortest wavelength) eg cosm(mx)cos(nz) showed energy increase.

[tomchor]

@glwagner Thanks for the detailed analyses and good points about the noise being grid-scale. That may indeed reduce the necessary time step for stability to unphysical values.

Also, just to be clear, I'm not saying that I think there's a bug in the code. I'm just saying that I, personally, wouldn't rule it out yet.

I guess part of my reasoning is due to the fact that I arrived at this MWE because of unphysical oscillations like these emerging due to physical (and as far as I can tell well-resolved) instabilities in my simulations. In these simulations the instabilities propagated into the stably-stratified, quasi-quiescent regions of the flow, and amplified similarly to the ones that the MWE above reproduces.

Of course the fact that these behaviors look alike to the naked eye doesn't prove they are indeed the same phenomenon and, like @glwagner suggested, the fact that we're starting with grid-scale noise here possibly matters.

If that's okay I'll try to investigate this a little more with grid-resolved noise before we move/close the issue. If the behavior goes away when everything is well-resolved, then I'm happy to close the issue.

[glwagner]

Rather than closing the issue I think it could be nice to turn this into a discussion, since it contains some advice about well-posed model setups with stable stratification (ie the need for small scale diffusivity).

[glwagner]

Also this doesn't blow up, right? I think with increasing energy, eventually nonlinearities and numerical diffusion would kick in. The integrated KE might exhibit oscillatory behavior on long time-scales. (The example in the OP only blows up because the TimeStepWizard increases the time step too much.)

[francispoulin]

Thanks everyone for the interesting discussion! I agree that this is of interest to many users and hope others benefit from it.

I am glad that @tomchor mentioned APE. Looking at the KE is very good but @glwagner has shown that the stratification is being modified, whichi suggests to me that APE might be important here. Are we able to look at KE + APE? I remember a discussion on this from a while ago and there was interest but I didn't follow up.

[tomchor]

I haven't looked at APE because it's harder to compute. Although maybe with the area and distance operators it wouldn't be too hard to calculate it via integral (Winters et al., 1995):

where $z*$ is the sorted $z$.

Note that initially the APE is zero since there are no fluctuations on top of the linear stable stratification. So the initial KE is indeed all the energy available to the system.

[glwagner]

I think the initial big drop has to do with energy being quickly exchange between KE and APE in this case. So I think that big drop is expected given the initial small-scale noise we include in the IC.

Perhaps but it may also be numerical diffusion. When the amplitude of the noise is increased, the drop is larger. When the amplitude of the initial noise is increased enough, the solution becomes smooth at long times and is stable. The relative amount of energy diffused is proportional to its amplitude because the numerical diffusivity is nonlinear (much like an LES diffusivity) in the amplitude of the solution.

[tomchor]

Okay so indeed I think this issue only happens due to poor resolution of the initial noise. I cooked up a MWE where I fixed the random seed, then used that to sprinkle 256 Gaussians (that are also horizontally periodic) throughout the domain. I did this because I wanted to have the "same noise", but change it from being under-resolved to well-resolved.

I ran this for a few cases ranging from poorly resolved (1 to 2 grid-points per gaussian) to well-resolved (about 8 grid-points per gaussian). I'm plotting three of the animations I got in order or increasing resolution (2, 4, and 8 grid-points per gaussian):

test_resnoise1.mp4

test_resnoise2.mp4

test_resnoise4.mp4

We can clearly see that the odd behavior decreases and then finally goes away as we start resolving the initial fluctuations more and more.

Some things caught my eye here though. Firstly, this is an example of the long-time solution depending on the initial conditions, which is really odd to me. Although I guess I shouldn't be so surprised, since this isn't physical turbulence, it's something else (and it's poorly resolved), but still I think it's interesting.

Also as the time progresses, even in the most resolved simulation, some grid-scale noise starts to emerge (which probably has to do with the advection scheme?). However, this grid-scale noise, contrary to the initial grid-scale noise, doesn't cause any trouble. I wonder why...

In any case, @glwagner please feel free to convert this issue to a discussion and thanks for the help!

[tomchor]

Also, here's the code I used for the simulations above, in case anyone's interested:

using Oceananigans
using Oceananigans.Units
using Random
Random.seed!(722)

grid = RectilinearGrid(topology=(Periodic, Flat, Bounded), size=(64, 32), x=(0, 200), z=(0, 100))

model = NonhydrostaticModel(; grid,
                            timestepper = :RungeKutta3,
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b)

# A kind of convoluted way to create x-periodic, resolved initial noise
σx = .5grid.Δxᶜᵃᵃ # x length scale of the noise
σz = .5grid.Δzᵃᵃᶜ # z length scale of the noise

N = 2^8 # How many Gaussians do we want sprinkled throughout the domain?
x₀ = grid.Lx * rand(N); z₀ = grid.Lz * rand(N) # Locations of the Gaussians

xₚ = x₀ .+ (grid.Lx .* [-2;;-1;;0;;1;;2]) # Make that noise periodic by "infinite" horizontal reflection
zₚ = reshape(repeat(z₀, 5), size(xₚ)) # Vertical direction isn't periodic, so we just need to repeat and reshape things 

resolved_noise(x, y, z) = sum(@. 1e-6 * exp(-(x-xₚ)^2/σx^2 -(z-zₚ)^2/σz^2))

constant_stratification(x, y, z) = 1e-5 * z
set!(model, u=resolved_noise, b=constant_stratification)

u, v, w = model.velocities

simulation = Simulation(model, Δt=20.0, stop_iteration=1e4)
progress(sim) = @info string("Iter: ", iteration(sim), ", time: ", prettytime(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

e = (u^2 + w^2) / 2
E = Field(Integral(e, dims=:))
outputs = merge(model.velocities, model.tracers, (; e, E))

filename = "test.jld2"
simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs;
                                                      filename,
                                                      schedule = IterationInterval(10),
                                                      overwrite_existing = true)

run!(simulation)

wt = FieldTimeSeries(filename, "w")
et = FieldTimeSeries(filename, "e")
bt = FieldTimeSeries(filename, "b")
Et = FieldTimeSeries(filename, "E")
times = wt.times
Nt = length(times)

using GLMakie
fig = Figure(resolution=(900, 450))

axw = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)", title="Vertical velocity")
axe = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)", title="Kinetic energy")
axb = Axis(fig[1, 3], xlabel="x (m)", ylabel="z (m)", title="Buoyancy")
axE = Axis(fig[2, 1:3], xlabel="Time (s)", ylabel="Volume averaged kinetic energy")

n = Observable(1)

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

x, y, z = nodes(et)

wlim = maximum(abs, wt)
elim = maximum(abs, et)

heatmap!(axw, x, z, wn, colorrange=(-wlim, wlim), colormap=:balance)
heatmap!(axe, x, z, en, colorrange=(0, elim), colormap=:solar)
heatmap!(axb, x, z, bn, colorrange=(0, 1e-3), colormap=:thermal)

t = @lift times[$n]
lines!(axE, times, Et[:])
vlines!(axE, t)

display(fig)

record(fig, "test.mp4", 1:Nt) do nn
    @info "Recording frame $nn of $Nt"
    n[] = nn
end