Velocity blowup in channel model?

[masonrogers14]

I'm trying to code a simple 2D channel model in v0.17.0 where the periodic dimension is the flat dimension. I keep experiencing seemingly spontaneous velocity blowup during spin-up despite changing parameters/removing most forcing. Here is an example:

Nx = 1
Ny = 256
Nz = 32
Lx = 1000e3
Ly = 1000e3
Lz = 4e3
N² = 1e-5
Δp = 1e-5 #will set to be a function of y
end_time = 1day

# ## Define a forcing function
pressure_gradient(x, y, z, t) = Δp / Lx
pressure_gradient_forcing = SimpleForcing(pressure_gradient)

## Instantiate the model
model = ChannelModel(
                   grid = RegularCartesianGrid(size = (Nx, Ny, Nz), x = (0, Lx), y = (-Ly/2, Ly/2), z=(-Lz, 0)),
                closure = ConstantIsotropicDiffusivity(ν=100, κ=100),
               coriolis = FPlane(f=1e-5),
                tracers = (:b,),
               buoyancy = BuoyancyTracer(),
                forcing = ModelForcing(v=pressure_gradient_forcing)
)

## Set initial condition.
b₀(x, y, z) = N² * z
set!(model, b=b₀)

## A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=90.0)

[ali-ramadhan]

Thanks for reporting this @masonrogers14!

In this case the pressure_gradient(x, y, z, t) = Δp / Lx forcing function may keep unrealistically accelerating u forever but I think the blowup is not related to that as we tried other forcing functions too. And decreasing the magnitude of Δp did not slow the blow up down.

I'm somewhat concerned that this issue may be related to PRs #572 and #573 so we should definitely investigate.

[glwagner]

Whats the stability condition associated with diffusivity for this problem?

[glwagner]

I get

julia> Δz = 4e3 / 32
125.0

julia> ν = 100
100

julia> Δt = Δz^2 / ν
156.25

So your time-step has to be substantially smaller than 156 s to be stable. 90 s may be pushing it... does reducing the maximum time-step size stabilize the problem?

[ali-ramadhan]

Yes we tried reducing the time step to a constant Δt=1 and Δt=1e-2 and it still blows up in the same number of iterations, suggesting to me that it wasn't a time-stepping issue.

It blows up in ~20 iterations so wizard.Δt is still pretty small.

[ali-ramadhan]

@glwagner suggested that recomputing w from continuity could be an issue here. We should try time-stepping w explicitly and running the example again.

[masonrogers14]

Reducing time steps for the example I shared solved the problem but fails for a more authentic example. Here is a better example, where I (very loosely) have a 2D y-z model with x-wind forcing in the southern hemisphere:

Nx = 1
Ny = 256
Nz = 32
Lx = 1000e3
Ly = 1000e3
Lz = 4e3
N² = 1e-5
end_time = 1day

## Wind stress BC
wind_stress(x,y,t) = -1e-5*sin(2π*y/Ly) * (y<0)
wind_stress_bf = BoundaryFunction{:z, Cell, Cell}(wind_stress)
velocity_bcs = ChannelBCs(top=BoundaryCondition(Flux, wind_stress_bf))

## Instantiate the model
model = ChannelModel(
                   grid = RegularCartesianGrid(size = (Nx, Ny, Nz), x = (0, Lx), y = (-Ly/2, Ly/2), z=(-Lz, 0)),
                closure = ConstantIsotropicDiffusivity(ν=100, κ=100),
               coriolis = BetaPlane(rotation_rate=1/day, latitude=0, radius=6400e3), #2π/day?
                tracers = (:b,),
               buoyancy = BuoyancyTracer(),
    boundary_conditions = BoundaryConditions(u=velocity_bcs)
)

## Set initial condition.
b₀(x, y, z) = N² * z
set!(model, b=b₀)

## A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=.001)

@glwagner -- explicitly timestepping w unfortunately did not solve the problem, but I'll keep looking into it

[glwagner]

Don't know if this is the problem, but the wind stress BoundaryFunction should be specified at Face, Cell; eg change the line with wind_stress_bf to

wind_stress_bf = BoundaryFunction{:z, Face, Cell}(wind_stress)

This is because the u velocity component is located at (Face, Cell, Cell) in (x, y, z) on the staggered grid.

[glwagner]

Note there are also a few changes to the script that would need to be made if the intention is to run on the GPU.

[masonrogers14]

Thanks @glwagner, I fixed the BoundaryFunction specification but am still getting blowup at iteration 23. I don't currently intend to run on a GPU but will keep that in mind

[masonrogers14]

I think I might be narrowing down the issue to implementation of boundary conditions. Perhaps I'm misunderstanding how ChannelModel works? The following code works for at least 10 days:

Nx = 1
Ny = 256
Nz = 32
Lx = 1000e3
Ly = 1000e3
Lz = 4e3
N² = 1e-5
end_time = 10day

# ## Creating boundary conditions
wind_stress(x,y,t) = -1e-5*sin(2π*y/Ly) * (y<0)
wind_stress_bf = BoundaryFunction{:z, Face, Cell}(wind_stress)
velocity_bcs = HorizontallyPeriodicBCs(top=BoundaryCondition(Flux,wind_stress_bf))
buoyancy_bcs = HorizontallyPeriodicBCs(bottom=BoundaryCondition(Gradient, N²))

## Instantiate the model
model = Model(
                   grid = RegularCartesianGrid(size = (Nx, Ny, Nz), x = (0, Lx), y = (-Ly/2, Ly/2), z=(-Lz, 0)),
                closure = ConstantIsotropicDiffusivity(ν=100, κ=100),
               coriolis = FPlane(f=1e-5),
                tracers = (:b,),
               buoyancy = BuoyancyTracer(),
    boundary_conditions = BoundaryConditions(b=buoyancy_bcs, u=velocity_bcs)
)

## Set initial condition.
b₀(x, y, z) = N² * z
set!(model, b=b₀)

## A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=20)

It fails after 10 iterations or so with .001s time steps when the Model is replaced with a ChannelModel and the HorizontallyPeriodicBCs are replaced with ChannelBCs.

[ali-ramadhan]

Holy shit this is a horrible alias!

https://github.com/climate-machine/Oceananigans.jl/blob/588890004e69cfc7db10472b12a9840b8a9ad7b6/src/boundary_conditions.jl#L339

[ali-ramadhan]

Yeah so you were correctly setting it up but apparently this alias meant that we were unknowingly applying horizontally periodic boundary conditions which would be incompatible with the pressure solver. Makes sense that it blew up very quickly no matter what.

Changing to

 boundary_conditions = ChannelSolutionBCs(b=buoyancy_bcs, u=velocity_bcs)

fixed it for me and now your example time steps fine on my laptop.

[ali-ramadhan]

This is the box model modified to use the correct ChannelSolutionBCs that I'm able to time step:

using Printf

using Oceananigans
using Oceananigans.Diagnostics

using Oceananigans: Cell, Face

Nx = 1
Ny = 256
Nz = 32
Lx = 1000e3
Ly = 1000e3
Lz = 4e3
N² = 1e-5
end_time = 10day

# ## Creating boundary conditions
wind_stress(x,y,t) = -1e-5*sin(2π*y/Ly) * (y<0)
wind_stress_bf = BoundaryFunction{:z, Face, Cell}(wind_stress)
velocity_bcs = ChannelBCs(top=BoundaryCondition(Flux,wind_stress_bf))
buoyancy_bcs = ChannelBCs(bottom=BoundaryCondition(Gradient, N²))

## Instantiate the model
model = Model(
                   grid = RegularCartesianGrid(size = (Nx, Ny, Nz), x = (0, Lx), y = (-Ly/2, Ly/2), z=(-Lz, 0)),
                closure = ConstantIsotropicDiffusivity(ν=100, κ=100),
               coriolis = FPlane(f=1e-5),
                tracers = (:b,),
               buoyancy = BuoyancyTracer(),
    boundary_conditions = ChannelSolutionBCs(b=buoyancy_bcs, u=velocity_bcs)
)

## Set initial condition.
b₀(x, y, z) = N² * z
set!(model, b=b₀)

## A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=20.0)

cfl = AdvectiveCFL(wizard)
dcfl = DiffusiveCFL(wizard)

while model.clock.time < end_time
    update_Δt!(wizard, model)

    time_step!(model; Δt=wizard.Δt, Nt=10)

    u, v, w = model.velocities
    u_max = maximum(abs, interior(u))
    v_max = maximum(abs, interior(v))
    w_max = maximum(abs, interior(w))

    t = model.clock.time
    @printf("Time: %s, Δt: %.2e CFL: %.2e, dCFL: %.2e, max (u, v, w): (%.2e, %.2e, %.2e)\n",
            prettytime(model.clock.time), wizard.Δt, cfl(model), dcfl(model), u_max, v_max, w_max)
end

[ali-ramadhan]

@masonrogers14 Just expanding on an earlier note by @glwagner on max_Δt: as diffusion is time stepped explicitly in Oceananigans, the maximum time step is limited by the stability criterion for the heat equation (see Von Neumann stability analysis for a derivation)

Δt < Δz^2 / 2κ

In practice you'll want Δt to be significantly below that limit, by a factor of 5x or more (adding each extra dimension to the problem lowers the limit by a factor of 2 I believe).

[masonrogers14]

@ali-ramadhan Thanks so much! I'm inclined to say the problem is resolved and will do more tests to make sure. Also the stability condition makes sense; I wasn't sure what the time stepper handled and what one has to address explicitly in the model

[ali-ramadhan]

Awesome! I'll close this issue then.

Also the stability condition makes sense; I wasn't sure what the time stepper handled and what one has to address explicitly in the model

The TimeStepWizard will make sure your CFL number (advective stability criterion) is below e.g. 0.2 but won't do anything about the diffusive CFL, although we could modify it to take into account both advective and diffusive stability criteria.