Evolving perturbations vs total fields in Nonhydrostatic model

[francispoulin]

@m3azevedo and me are doing a simulation of an unstable baroclinic jet in the context of the Nonhydrostatic model and having some success but some difficulties as well.

We have tried two different approaches:

Method 1: define the jet and stratification to be part of the background field and evolve the perturbations.

Method 2: define the jet and stratification as part of the initial conditons and evolve the total field.

These are mathematically equivalent and should both yield the same results, but we are finding a difference and believe there is a problem with method 2 (evolving the total fields).

When pick there to be zero perturbations, the jet is an exact solution to the governing equations (even though it is an unstable one) and it should persist, until perturbations arise because of numerical errors.

The good news is that when we use method 1, evolve the perturbations, when we pick zero perturbations, no pertrubations develop for a long time.

The not so good news, is that when we use method 2, evolve the total field, a distrubance develops right away,in the center of the jet that spans the entire depth, and this seems to radiate gravity waves to the left and right. The amplitude of the distrubance is small and it does decrease with increasing resolution but it seems problematic.

I am going to include an animation that shows the problem that develops in the first 6 hours and the two codes, which are almost identical. Is there something that we are doing wrong here? @glwagner ?

Method 1: Evolve perturbations

using Oceananigans
using Oceananigans.Units
using NCDatasets, Printf
using Statistics, Random
using LinearAlgebra: norm

## background consts
f     = 0.864e-4 
fₕ    = 0.0 
N²    = (3.7e-3)^2 
ν     = 0.36/2.2e5
Umax  = 0.36
Lx     = 10000
Ly     = 30000
Lⱼ    = 2000
Lz    = 1000
D     = 200
z0    = -Lz/2

grid = RectilinearGrid(
    CPU(),
    size=(1, 200, 100),
    x= (-Lx, Lx),
    y = (-Ly, Ly),
    z = (-Lz, 0),
    topology=(Periodic, Bounded, Bounded)
    )

U(x, y, z, t) = Umax / cosh(y/Lⱼ)^2 * exp(-(z-z0)^2/D^2)
B(x, y, z, t) = N² * z + 2*f*Umax*Lⱼ/D^2 * (tanh(y/Lⱼ)) * (z-z0) * exp(-(z-z0)^2/D^2)

Random.seed!(25)
uᵢ(x, y, z) = 0.0
bᵢ(x, y, z) = 0e-6*rand()

model = NonhydrostaticModel(; grid,
    background_fields = (u=U, b=B),
    coriolis = FPlane(f = f),
    buoyancy = BuoyancyTracer(),
    tracers  = (:b,),
    advection = WENO(grid=grid),
    closure = VerticalScalarDiffusivity(ν=ν, κ=ν)
)
set!(model, u = uᵢ, b = bᵢ)

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

function progress(sim)
    umax = maximum(abs, sim.model.velocities.u)
    bmax = maximum(abs, sim.model.tracers.b)
    @info @sprintf("Iter: %d, time: %.2e, max|u|: %.2e, max|b|: %.2e",
       iteration(sim), time(sim), umax, bmax)

    return nothing
end
simulation = Simulation(model; Δt=30, stop_time=6hours)
simulation.callbacks[:p] = Callback(progress, TimeInterval(1minutes))

u, v, w = model.velocities
b = model.tracers.b
outputs = (; v, u, w, b)

simulation.output_writers[:fields] = NetCDFOutputWriter(
        model, outputs;
        filename = "NH_BC_jet_pert_fields.nc",
        schedule = TimeInterval(1minutes),
        array_type = Array{Float32},
        overwrite_existing = true)

perturbation_norm(args...) = norm(u)

simulation.output_writers[:growth] = NetCDFOutputWriter(
        model, (; perturbation_norm),
        filename = "NH_BC_jet_pert_norm.nc",
        schedule = TimeInterval(1minutes),
        dimensions = (; perturbation_norm = ()),
        overwrite_existing = true)

run!(simulation)

Method 2: evolve total field

using Oceananigans
using Oceananigans.Units
using NCDatasets, Printf
using Statistics, Random
using LinearAlgebra: norm

## background consts
f     = 0.864e-4 
fₕ    = 0.0 
N²    = (3.7e-3)^2 
ν     = 0.36/2.2e5
Umax  = 0.36
Lx     = 10000
Ly     = 30000
Lⱼ    = 2000
Lz    = 1000
D     = 200
z0    = -Lz/2

grid = RectilinearGrid(
    CPU(),
    size=(1, 200, 100),
    x= (-Lx, Lx),
    y = (-Ly, Ly),
    z = (-Lz, 0),
    topology=(Periodic, Bounded, Bounded)
    )

U(x, y, z) = Umax / cosh(y/Lⱼ)^2 * exp(-(z-z0)^2/D^2)
B(x, y, z) = N² * z + 2*f*Umax*Lⱼ/D^2 * (tanh(y/Lⱼ)) * (z-z0) * exp(-(z-z0)^2/D^2)

Random.seed!(25)
uᵢ(x, y, z) = U(x, y, z)
bᵢ(x, y, z) = B(x, y, z) + 0e-6*rand()

b_bc = FieldBoundaryConditions(top = GradientBoundaryCondition(N²), bottom = GradientBoundaryCondition(N²))

model = NonhydrostaticModel(; grid,
    coriolis = FPlane(f = f),
    buoyancy = BuoyancyTracer(),
    tracers  = (:b,),
    advection = WENO(grid=grid),
    closure = VerticalScalarDiffusivity(ν=ν, κ=ν),
    boundary_conditions = (b = b_bc, )
)
set!(model, u = uᵢ, b = bᵢ)

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

function progress(sim)
    umax = maximum(abs, sim.model.velocities.u)
    bmax = maximum(abs, sim.model.tracers.b)
    @info @sprintf("Iter: %d, time: %.2e, max|u|: %.2e, max|b|: %.2e",
       iteration(sim), time(sim), umax, bmax)

    return nothing
end
simulation = Simulation(model; Δt=30, stop_time=1hours)
simulation.callbacks[:p] = Callback(progress, TimeInterval(1minutes))

u, v, w = model.velocities
b = model.tracers.b
outputs = (; v, u, w, b)

simulation.output_writers[:fields] = NetCDFOutputWriter(
        model, outputs;
        filename = "NH_BC_jet_fields.nc",
        schedule = TimeInterval(1minutes),
        array_type = Array{Float32},
        overwrite_existing = true)

perturbation_norm(args...) = norm(u)

simulation.output_writers[:growth] = NetCDFOutputWriter(
        model, (; perturbation_norm),
        filename = "NH_BC_jet_norm.nc",
        schedule = TimeInterval(1minutes),
        dimensions = (; perturbation_norm = ()),
        overwrite_existing = true)

run!(simulation)

[francispoulin]

https://github.com/CliMA/Oceananigans.jl/assets/8239041/3e75ec65-3a2b-4063-acf0-423b016d63fe

[simone-silvestri]

I probably would have to think a bit about it, but aren't baroclinic instabilities developing on the x-y plane? Can you represent a baroclinic instability process with 1 point in the jet direction?

[francispoulin]

This is a baroclinic jet and you are correct that having a Flat topology prevents any baroclinic instabilities from happening. We are actually studying inertial instabilities, and they can occur in a y-z plane.

[francispoulin]

One thing that might be easier to ask for help with his, I wanted to print out the norm of the perturbation buoyancy in both cases. When we evolve the perturbation this is easy, but when we evolve the total field, I couldn't figure out how to do it. Do you have a recommendation?

This would be helpful to determine how the norm changes with resolution.

[simone-silvestri]

I gues something like this would work

B₀ = CenterField(grid)
set!(B₀, B)
b = model.tracers.b
b′² = Field(Average((b - B₀)^2))

[francispoulin]

Thanks for the suggestion @simone-silvestri .

Building on your suggestion, what I tried below seemd to do the trick.

Cheers!

B_background = CenterField(grid)
set!(B_background, B)

function progress(sim)
    umax = maximum(abs, sim.model.velocities.u)
    bmax = maximum(abs, sim.model.tracers.b - B_background)
    @info @sprintf("Iter: %d, time: %.2e, max|u|: %.2e, max|b|: %.2e",
       iteration(sim), time(sim), umax, bmax)

    return nothing
end

[francispoulin]

@simone-silvestri , thanks to your help I am now able to print out the norms of the perturbation buoyancy in the case when we evolve the total fields.

The error that arises in this case, but not when we evolve the perturbation fields, decreases quadratically with increasing resolution. I could pick a resolution where this is small but in the first hour the norm or buoyancy grows by three orders of magnitude, which looks like an instability. It is not noisy and very well resolved, it can't be a baroclinic instability and does not look anything like an inertial instability.

[glwagner]

Is it because your initial condition is only geostrophically balanced in the continuous limit? You may need to compute the velocity field that's in discrete geostrophic balance with your initial buoyancy field.

You are using analytical formula for both U and B:

U(x, y, z) = Umax / cosh(y/Lⱼ)^2 * exp(-(z-z0)^2/D^2)
B(x, y, z) = N² * z + 2*f*Umax*Lⱼ/D^2 * (tanh(y/Lⱼ)) * (z-z0) * exp(-(z-z0)^2/D^2)

but to satisfy the discrete geostrophic balance you would have to compute one of those fields from the other one, in a manner consistent with how Coriolis forces are implemented.

[glwagner]

Something like this might work:

B(x, y, z) = N² * z + 2*f*Umax*Lⱼ/D^2 * (tanh(y/Lⱼ)) * (z-z0) * exp(-(z-z0)^2/D^2)
set!(model, b=B) # this calls `update_state!`, and thus calculates the hydrostatic pressure associated with buoyancy

f = model.coriolis.f
ph = model.pressures.pHY′
U =  - ∂y(ph) / f
set!(model, u=U)

Curious to see if that solves it.

And just for sanity we can look into the source to ensure that the hydrostatic pressure is indeed calculated during set!:

https://github.com/CliMA/Oceananigans.jl/blob/4a31c0cd15e4d3e414decde7ba1bb285a1d015f0/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl#L49-L54

Not a terrible idea to plot model.tracers.b, ph, and U too (for that you'll need to write U = compute!(- ∂y(ph) / f)).

[francispoulin]

Thanks @glwagner for the suggestion. We will try this out and see if we can improve our results, and if this instability disappears.

Actually, we found that the Hydrostatic model suffers from the same problem, and since it solves for the total field that is maybe not surprising. So confirming the hydrostatic pressure is a very good idea!

[francispoulin]

@glwagner : I was able to get your suggestion working very easily, thanks again.

Unfortuantely it doesn't help is removing the error that is introduced in the centre or the fact that it grows exponentially (I am guessing).

I have looked at the pressure and it looks smooth but I suppose looking at the pertrubation pressure might be more telling.

[glwagner]

Can you confirm that there is hydrostatic balance at the discrete level? It might help to plot the deviation from thermal wind, for example, to get a clue as to what's going on.

[glwagner]

PS can we convert this to a discussion? I very much doubt this is a bug in the source code and it might be nice to have the discussion format to organize our ideas to solve this problem

[navidcy]

PS can we convert this to a discussion? I very much doubt this is a bug in the source code and it might be nice to have the discussion format to organize our ideas to solve this problem

I was thinking the same!

[m3azevedo]

Thank you for the great suggestions so far. I see no problem with converting this thread to a discussion.

Although the problem is not solved yet, we have done some investigation around it: 1) Increasing the precision of the output to Float64 eliminates the initial error seen at the first few frames of the animation.

2) Using a barotropic version of the jet

U(x, y, z) = Umax / cosh(y/Lⱼ)^2
B(x, y, z) = N² * z

leads to no perturbation when evolving the total fields.

3) Using a discretized version of the thermal wind balance is a sound suggestion, but it does not explain why the error only appears when evolving the total fields.

[glwagner]

Using a discretized version of the thermal wind balance is a sound suggestion, but it does not explain why the error only appears when evolving the total fields.

Why not? When evolving the total fields, the condition of geostrophic balance is computed during the prognostic evolution of the total variables. When evolving just the perturbation, the algorithm for mean-perturbaiton decomposition assumes that the mean fields are balanced. There is no need to worry about the tiny difference between discrete and continuous balance in the case that you only evolve perturbations. But if you want to evolve the geostrophically balanced state, you must ensure that it is perfectly balanced to within machine precision, or you will have an adjustment. That's the basis for my suggestion. I'm not sure the computation I recommended achieves numerical balance (something might be missing) --- that should be checked.

[glwagner]

Increasing the precision of the output to Float64 eliminates the initial error seen at the first few frames of the animation.

I'm not sure I understand. The error is only associated with post-processing and not with the actual simulation?

[francispoulin]

@glwagner : did you want to convert this to a discussion before we go further in our chat?

[glwagner]

Yes if youre ok with that I will