Wind-forced simulation is exploding with lateral noise

[iuryt]

(UPDATE: it is worth reading all @simone-silvestri explanations, but basically, I was dumb enough to integrate N2 in the wrong way and initialize with unstable stratification! haha)

I am running an experiment to simulate a classical internal wave generation problem and need your help.

The idea is to start with the ocean at rest in an equatorial-beta simulation that is solved only in (y, z) and add a sine wind stress in a region 2750km north of the equator. The wind time-variation is Gaussian with a width of three standard deviations of two inertial periods. The idea is to leave unbalanced currents that will oscillate inertially.

In the video below I show the horizontal velocities (u on the left and v on the right) in an approximate region of the storm. The top panels are horizontal profiles of u and v at 20m and the right panel shows vertical profiles of u and v in place of maximum stress.
I apologize for the preliminary plots with labels covering each other and no time tag, but the whole video is for 10h of simulation time.

For some reason the simulation is blowing up due to horizontal oscillations (please note that the vertical velocity profiles are smooth, but the horizontal profiles at 20 m in the upper panels are quite noisy, especially for v).

I'm under the impression that it's something very obvious, but I swear I tried everything I thought I could improve and I couldn't. I just can't do it. So I'm posting here to see if anyone has an idea.

What I tried so far:

• change horizontal resolution: now using 1km but changed from 200m to 2km.
• change vertical resolution. Now using 100 points, changed to 1000 points. still blow up
• change biharmonic diffusivity from 1e3 to 1e8. A little bit better for higher diffusivity, but blow up sooner for 1e8
• decreasing maximum depth to 1000 gets a little bit better, but still blow up.
• I also tried to use biharmonic horizontal divergence damping and still blow up

I also think that it is mixing too much in the vertical as buoyancy lines seem to go up, but I don't think this is the reason why this is blowing up.

Any suggestions?

By the way, what is the best way to also save the boundary conditions in the output file? I wanted to save the wind stress.

movie_custom.mp4

using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity, HorizontalDivergenceFormulation
using Oceananigans.Advection: VelocityStencil

#--------------- Grid

const Ny = 6000 # number of points in y
const Nz = 100 # number of points in z
const H = 5000meters # maximum depth
const yc = 2000

grid = RectilinearGrid(GPU(),
    size=(Ny, Nz),
    halo=(3,3),
    y=((yc-3000)kilometers, (yc+3000)kilometers), 
    z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters,
    topology=(Flat, Bounded, Bounded)
)

#--------------- Coriolis

coriolis = BetaPlane(β=2.3e-11, latitude=0) # equatorial beta plane

#--------------- Bottom drag

const cᴰ = 1e-4 # quadratic drag coefficient
const c = 2750kilometers     # center of the storm
const L = pi / (250kilometers) # horizontal wavenumber of the storm
const Ω = 7.27E-5
const latitude = c/112kilometers
const ip = 2pi/(2Ω*sin(pi*latitude/180))

# const u₁₀ = 1    # m s⁻¹, average wind velocity 10 meters above the ocean
# const csᴰ = 2.5e-3 # dimensionless drag coefficient
# const ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
# const ρₒ = 1025 # kg m⁻³, average density of water
# @inline U(x, y, t) = ifelse(abs(L*(y-c))<pi, u₁₀*sin(L*(y-c))*exp(-((t-1.5*ip)^2)/(2*(2*ip/3)^2))*exp(-(y-c)^2/(2*(pi/L/3)^2)), 0)
# @inline Qᵘ(x, y, t) = - ρₐ / ρₒ * csᴰ * U(x, y, t) * abs(U(x, y, t)) # m² s⁻²
@inline Qᵘ(x, y, t) = 1e-4*sin(L*(y-c))*exp(-((t-1.5*ip)^2)/(2*(2*ip/3)^2))*exp(-(y-c)^2/(2*(pi/L/4)^2))

surface_drag_bc_u = FluxBoundaryCondition(Qᵘ)

@inline bottom_drag_u(x, y, t, u, v, cᴰ) = - cᴰ * u * sqrt(u^2 + v^2)
@inline bottom_drag_v(x, y, t, u, v, cᴰ) = - cᴰ * v * sqrt(u^2 + v^2)

bottom_drag_bc_u = FluxBoundaryCondition(bottom_drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
bottom_drag_bc_v = FluxBoundaryCondition(bottom_drag_v, field_dependencies=(:u, :v), parameters=cᴰ)

u_bcs = FieldBoundaryConditions(bottom = bottom_drag_bc_u, top = surface_drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = bottom_drag_bc_v)


#--------------- Sponges

const sponge_size = 500kilometers
const slope = 40kilometers
const ymin = minimum(ynodes(Center, grid))
const ymax = maximum(ynodes(Center, grid))

@inline mask_func(x,y,z) = ((
     tanh((y-(ymax-sponge_size))/slope)
    *tanh((y-(ymin+sponge_size))/slope)
)+1)/2

mom_sponge = Relaxation(rate=1/30minutes, mask=mask_func, target=0)


#--------------- Closure

boundary_layer_turbulence_closure = RiBasedVerticalDiffusivity()

# horizontal_viscosity = ScalarBiharmonicDiffusivity(HorizontalDivergenceFormulation(), ν=1e6, κ=1e6)
horizontal_viscosity = HorizontalScalarBiharmonicDiffusivity(ν=1e6, κ=1e6)

# vertical_viscosity   = VerticalScalarDiffusivity(ν=1e-5, κ=1e-5)

closures = (horizontal_viscosity, boundary_layer_turbulence_closure)

#--------------- Instantiate Model

model = NonhydrostaticModel(grid = grid,
                            advection = WENO5(vector_invariant = VelocityStencil()),
                            coriolis = coriolis,
                            closure = closures,
                            forcing = (u=mom_sponge, v=mom_sponge, w=mom_sponge),
                            tracers = (:b), 
                            buoyancy = BuoyancyTracer(),
                            boundary_conditions = (u=u_bcs, v=v_bcs)
)

#--------------- Initial Conditions

const a = 1e-5 # stratification at the pycnocline
const b = 0    # stratification below pycnocline
const h = 50   # mixed layer depth
const ht = 30  # height of the transition layer above pycnocline
const hp = 500 # height of the pycnocline

@inline N2(x,y,z) = ifelse( z >= -h-ht, a*(exp(-(z+h+ht)^2/(2*(ht/3)^2))), ((a-b)*(exp(-(z+h+ht)^2/(2*(hp/3)^2)))+b))

N2_field = Field{Center, Center, Center}(grid)
set!(N2_field, N2)

# include function for cumulative integration
include("cumulative_vertical_integration.jl")

b_field = 9.82 + cumulative_vertical_integration!(N2_field)

set!(model; b = b_field)

#--------------- Simulation and outputs

simulation = Simulation(model, Δt = 5minutes, stop_time = 30hours)

wizard = TimeStepWizard(cfl=1.0, max_change=3, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

ui,vi,wi = model.velocities
bi = model.tracers.b

outputs = (;
    b=@at((Center, Center, Center), bi),
    u=@at((Center, Center, Center), ui),
    v=@at((Center, Center, Center), vi),
    w=@at((Center, Center, Center), wi),
    N2=@at((Center, Center, Center), ∂z(bi)),
    Ri=@at((Center, Center, Center), ∂z(bi)/(∂z(ui)^2+∂z(vi)^2)),
)


simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, outputs, filename = "../data/output_wind.nc",
                        overwrite_existing=true,
                        schedule=TimeInterval(10minutes))

#--------------- Printing Progress

using Printf

function print_progress(simulation)
    u, v, w = simulation.model.velocities

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(print_progress, TimeInterval(10minutes))

run!(simulation)

The file cumulative_vertical_integration.jl

using KernelAbstractions: @index, @kernel
using KernelAbstractions.Extras.LoopInfo: @unroll
using Oceananigans.Architectures: device_event, architecture
using Oceananigans.Utils: launch!
using Oceananigans.Grids
using Oceananigans.Operators: Δzᶜᶜᶜ

@kernel function _cumulative_vertical_integration!(field, grid)
    i, j = @index(Global, NTuple)

    integ = 0
    @unroll for k in grid.Nz : -1 : 1
        integ = integ + @inbounds field[i, j, k] * Δzᶜᶜᶜ(i, j, k, grid)
        @inbounds field[i, j, k] = integ
    end
end

function cumulative_vertical_integration!(input_field)
    field = Field(input_field)
    grid = field.grid
    arch = architecture(grid)

    event = launch!(arch, grid, :xy,
                    _cumulative_vertical_integration!, field, grid,
                    dependencies = device_event(arch))

    wait(device_event(arch), event)

    return field
end

[simone-silvestri]

Two comments

1. vector invariant advection is good for 3D simulations but not so much for 2D because it does not implicitly diffuse the vertical direction. For Nohydrostatic simulations, vector invariant mode is not supported, so WENO(vector_invariant = VelocityStencil()) falls back to WENO()
2. Try using a simple vertical diffusion parametrization VerticalScalarDiffusivity(ν = 1e-4, κ = 1e-5) + ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0)

Be careful with the advective and diffusive CFL. For biharmonic diffusion ν < Δ⁴ / 32Δt, which means that if you are time stepping with Δt = 5minutes, you are limited to ν < 1e8. Having such a large biharmonic diffusion might be the whole issue here.
You also see that your simulation crashes when you reach u ≈ v ≈ 1e-3. If you check your w velocity, it might be that you are violating the vertical CFL condition. Also, you are using an Adams Bashfort time-stepping method which is quite restrictive in terms of time step

I would just try for the moment a simulation with no closure, a runge kutta time stepping method (timestepper = :RungeKutta3 which allows larger CFL than the default Adams Bashforth) and a CFL of 0.75

let us know how it goes!

[iuryt]

My God! How did I make this stupid mistake? haha I remember checking it before, but I guess because I'm used to density profiles I thought it was all right. Sorry for taking your time with something so simple that I should have checked it more carefully myself. Now everything is working fine.

[iuryt]

@simone-silvestri
Here is an updated video for the problem. But I have a question. How can we easily save the boundary conditions as an output?

movie_custom.19.mp4

[simone-silvestri]

Hmm, this looks indeed a bit noisy. Maybe a biharmonic viscosity would help.
To save the BC, you could define a field that contains it.

BC = Field((Face, Center, Nothing), grid)

outputs = (;
    b=@at((Center, Center, Center), bi),
    u=@at((Center, Center, Center), ui),
    v=@at((Center, Center, Center), vi),
    w=@at((Center, Center, Center), wi),
    N2=@at((Center, Center, Center), ∂z(bi)),
    Ri=@at((Center, Center, Center), ∂z(bi)/(∂z(ui)^2+∂z(vi)^2)),
    BC
)

Then you can calculate it in a callback

function update_bc(simulation)
  @inline Qᵘtemp(x, y) = Qᵘtemp(x, y, simulation.model.clock.time)
  set!(BC,  Qᵘtemp)
end

simulation.callbacks[:update_bc] = Callback(update_bc, TimeInterval(10minutes))

[iuryt]

Thanks!

[glwagner]

You could probably also use KernelFunctionOperation to avoid using the callback