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]

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

[glwagner]

I'd be curious what happens when other numerical methods (like spectral methods) simulate grid scale inviscid stratified flows.

[francispoulin]

I have used spectral methods to simulate nearly inviscid flow and can say you always need a fitler or hyperviscosity to control the build up at small scales. There is no closure in this problem so maybe that is part of the problem? Though I would have thought that the 5th order upwinding would have enough dissipation to remove the grid scale noise.

This example does seem to suggest that using random initial conditions without any closure could be problematic. Do people agree? Maybe I am missing something.

[glwagner]

Here, upwinded advection has little effect on very small amplitude solutions that are dominated by linear dynamics. I guess filtering for spectral methods removes a fixed fraction of energy at a given wavenumber so doesn't require nonlinear dynamics to act. For sure adding some explicit diffusion is the typical way to "filter" energy in a finite volume context. Just curious whether errors in psuedospectral method would also grow in a similar way (ie increasing KE...)

[francispoulin]

Interesting. Spectral methods don't differentiate between linear and nonlinear dynamics, so that is one difference. They just filter the small scales. I have certainly seen a number of spectral methods have small scale features grow very large in amplitude, usually because the filtering wasn't strong enough. I don't remember whether it compared to this scenario but when spectral methods have small scales grow it's usually said that the filteirng isn't done correctly.