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]

b_yz_NH.mp4

[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!:

Oceananigans.jl/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl

Lines 49 to 54 in 4a31c0c

	if enforce_incompressibility
	FT = eltype(model.grid)
	calculate_pressure_correction!(model, one(FT))
	pressure_correct_velocities!(model, one(FT))
	update_state!(model)
	end

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

[francispoulin]

@glwagner : We have followed up on your suggestion and have computed the error in geostrophic and hydrostatic balance at the initial time. The animation included here shows the results for the first 10 minutes.

In the geostrophic equation we get errors on the order of O(1e-7), which is about 100 times smaller, than the two terms in the equation. This shows that the error is about 1 percent.

The hydrostatic equation, on the other hand, looks to be zero from the animation O(1e-17) but upon further investigation, we find that the pressure at the bottom is O(1e-2) and at the top we have O(1e-6). The pressure in the interior is zero to machine precision. Something has gone wrong in the calculation for the hydrostatic pressure at the very top and the very bottom, and wonder how we can fix this.

gbhb.mp4

[francispoulin]

I have thougth about this a bit more and I believe I have identified two sources of error.

First, when I look at the vertical derivative of pressure at the top and bottom it is 0. We are not imposing what boundary conditions to impose on pressure and I guess it's assuming homogeneous Neumann boundary conditions, which is not correct.

Is there something else we need to do after ∂z(p) to enforce correct boundary conditions for pressure?

Second, we are not being careful enough with the positioning of variables. To test geostrophic and hydrostatic balance we used the following equations,

geoeqn = ∂y(p) + f*u
hydeqn = ∂z(p) - b

Since both u and p are at centres in the y direction, computing a y derivative then moves pressure to a different location.

Similary, pressure and buoyancy are both at centres in the vertical, and the z derivative then moves one of the fields to the wrong location.

I presume we need to do an 2 dimensional interpolation after we compute the derivative in each case to end up at the centres?

@glwagner @simone-silvestri ?

[glwagner]

Hmm yes it is a challenge to compute u given pressure I realize now. There are a few discretizations of Coriolis. FPlane is here:

Oceananigans.jl/src/Coriolis/f_plane.jl

Line 45 in b99b467

@inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = coriolis.f * ℑxyᶜᶠᵃ(i, j, k, grid, U[1])

involving 2D interpolation in xy.

In general do we have to solve an elliptic equation?

[glwagner]

Buoyancy is relatively easier; we interpolate in z:

Oceananigans.jl/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl

Line 16 in b99b467

@inbounds pHY′[i, j, k] = pHY′[i, j, k+1] - z_dot_g_bᶜᶜᶠ(i, j, k+1, grid, buoyancy, C) * Δzᶜᶜᶠ(i, j, k+1, grid)

[francispoulin]

Thanks for thinking about this some more and the links.

No, I don't believe that we need to solve an elliptic equation in this context but maybe someone else things differently?

It occurs to me that we have analytical expressions for pressure, velocity and buoyancy. I think it would be better to use the analytical expression for pressure and compute the derivatives and interpolations directly to find the velocity and buoyancy. This should give us that the geostrophic and hydrostatic equations are exactly zero for any grid. Thermal wind might not be zero but would be happy to be wrong.

What do you think about this @glwagner? I will give it a try tomorrow but I know the interpolations can be a bit tricky.

[glwagner]

I think it would be better to use the analytical expression for pressure and compute the derivatives and interpolations directly to find the velocity and buoyancy.

But you can't compute either the velocity or buoyancy directly from hydrostatic pressure. You can define a pressure field and take the y-derivative, but in the discrete equations, this must be balanced by the interpolated velocity field. You thus don't have a discrete algebraic relation for velocity; the relation is non-local due to the interpolation (you could think of this as an integral equation).

The discrete balance relations are

f * Ixy(u) = dy(p)
- Iz(b) = dz(p)

where Ixy(u) is a non-local operator that interpolates u from u-locations to v locations (where dy(p) is located).

If you want to avoid those complications, I think maybe you can just try to increase the resolution higher and higher. I believe the out-of-balance component of the initial condition associated with replacing discrete expressions by continuous ones should get smaller and smaller...