Wind-driven Ocean Circulation

[iuryt]

Hi,

I am not sure this is my actual problem, but I am trying to run a very simple wind-driven currents model for the South Atlantic with the code here below.

using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity

const H = 1000 # maximum depth
const Nz = 48

underlying_grid = RectilinearGrid(GPU(), 
    size=(800, 400, Nz),
    x=(-4000kilometers, 4000kilometers),
    y=(-2000kilometers, 2000kilometers),
    z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters,
    topology=(Bounded, Bounded, Bounded),
    halo=(4, 4, 4)
)


@inline bathymetry(x,y) = (-1e3*(tanh((x-1e6*sin(pi*(y-0.5e6)/4e6)+2.5e6)/1e5)+1)/2)-2
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))


@inline u₁₀(x,y) = -10*sin(pi*y/2000kilometers)  # m s⁻¹, average wind velocity 10 meters above the ocean
const cᴰ = 2.5e-3 # dimensionless drag coefficient
const ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
const ρₒ = 1025.0 # kg m⁻³, average density at the surface of the world ocean

@inline Q(x,y,t) = - ρₐ / ρₒ * cᴰ * u₁₀(x,y) * abs(u₁₀(x,y)) # m² s⁻²


no_slip_bc = ValueBoundaryCondition(0.0)
u_bcs = FieldBoundaryConditions(no_slip_bc, top = FluxBoundaryCondition(Q))
v_bcs = FieldBoundaryConditions(no_slip_bc, top = FluxBoundaryCondition(nothing))


horizontal_closure = HorizontalScalarDiffusivity(ν=10, κ=10)
# vertical_closure = ScalarDiffusivity(ν=1e-5, κ=1e-5)
vertical_closure = RiBasedVerticalDiffusivity()

model = NonhydrostaticModel(; grid,
                            coriolis = BetaPlane(latitude = -35),
                            closure = (horizontal_closure,vertical_closure),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b,
                            advection = WENO(),
                            boundary_conditions = (u=u_bcs, v=v_bcs))

# starting with a constant buoyancy
set!(model, b=-9.8)


simulation = Simulation(model, Δt = 1minutes, stop_time = 10day)
wizard = TimeStepWizard(cfl=1.0, max_change=0.5, max_Δt=3hour)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))


b, = simulation.model.tracers
u, v, w = simulation.model.velocities

extra_outputs = (; 
    u=@at((Center, Center, Center), u),
    v=@at((Center, Center, Center), v),
    w=@at((Center, Center, Center), w),
    Zeta=@at((Center, Center, Center), ∂x(v)-∂y(u)),
)

filename = "../data/output_box.nc"

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, merge(model.tracers, extra_outputs),
                     overwrite_existing=true,
                     filename = filename,
                     schedule=TimeInterval(3hours))



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(1hour))


run!(simulation)

You guys can ignore the part of the code where I define the output and progress messages, but would like to have a general opinion on that. What happen is the the following code runs for a nonbathymetry case, but velocities goes NaN for this case.. If I put the shallowest level in -200m instead of 0, the code runs a little bit more.
Not sure what I am doing wrong here, but I suspect it is related to the immersed boundary and the boundary condition at the surface.

I am trying to create a zeroth-order simulation to specify later to a more real bathymetry and winds.
Right now I am running for homogeneous water column and simple sine winds.

Does anything pop up just by looking at the code?
Any suggestions?

[glwagner]

A few questions and comments:

• What does the bathymetry look like? Can you plot it? (Also, I suggest splitting up the bathymetry definition into more lines, because it is hard to understand what (-1e3*(tanh((x-1e6*sin(pi*(y-0.5e6)/4e6)+2.5e6)/1e5)+1)/2)-2 does at a glance)
• What vertical grid spacing is produced by z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters ?
• Can you try the same code with HydrostaticFreeSurfaceModel? You should definitely use the hydrostatic model for this case, which will be more accurate (the nonhydrostatic pressure solver is not exact in the presence of bathymetry)
• I suggest starting with the default free slip boundary conditions. Hard to say whether the no-slip condition being used in conjunction with RiBasedVerticalDiffusivity causes a problem. A more appropriate model for bottom drag would be quadratic bottom drag via FluxBoundaryCondition, and that also might be better behaved with RiBasedVerticalDiffusivity.
• (Side comment: I suggest imposing the momentum flux directly, which will reduce the number of parameters your script depends on and thus will help with reproducibility)

[simone-silvestri]

Are you sure you are using the same version of Julia? That line has been recently changed because we switched to Julia 1.9. which apparently defines fkeep! in another way than Julia 1.8, @navidcy may know more about it

[iuryt]

The newest version of Julia from the cluster I am using is 1.8.3. I will check with them if they can install Julia 1.9 and let you guys know.

[glwagner]

Ok, so that means that Oceananigans does not work for versions 1.8 and earlier, correct? We need to define a compat then. Right @navidcy ?

[iuryt]

Alright, so now I am using Julia 1.9.2 and Oceananigans 0.84.1

I am running the code below for only 10 minutes of simulation and it is generating this boundary effect. I am using the hydrostatic model and with the underlying grid (no bathymetry), and no bottom drag. I am setting up FluxBoundaryCondition(nothing) for free-slip conditions. Here I am plotting at yC=0 and zC=0 the velocities v (blue) and u (orange). Is this expected? Any suggestion? Maybe I should add a sponge?

The code for the simulation

using Oceananigans
using Oceananigans.Units
using Oceananigans.BoundaryConditions: ImpenetrableBoundaryCondition
using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity

const H = 1000 # maximum depth
const Nz = 48

underlying_grid = RectilinearGrid(GPU(), 
    size=(800, 400, Nz),
    x=(-4000kilometers, 4000kilometers),
    y=(-2000kilometers, 2000kilometers),
    z=(H * cos.(LinRange(π/2,0,Nz+1)) .- H)meters,
    topology=(Bounded, Bounded, Bounded),
    halo=(4, 4, 4)
)


#@inline bathymetry(x,y) = (-1e3*(tanh((x-1e6*sin(pi*(y-0.5e6)/4e6)+2.5e6)/1e5)+1)/2)-2
#grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))


const cbᴰ = 1e-4 # quadratic bottom drag coefficient
const cᴰ = 2.5e-3 # dimensionless surface drag coefficient
const ρₐ = 1.225  # kg m⁻³, average density of air at sea-level
const ρₒ = 1025.0 # kg m⁻³, average density at the surface of the world ocean

@inline u₁₀(x,y) = -10 * sin(pi * y / 2000kilometers)  # m s⁻¹, average wind velocity 10 meters above the ocean
@inline Q(x,y,t) = - ρₐ / ρₒ * cᴰ * u₁₀(x,y) * abs(u₁₀(x,y)) # m² s⁻²

#@inline bottom_drag_u(x, y, t, u, v, cbᴰ) = - cbᴰ * u * sqrt(u^2 + v^2)
#@inline bottom_drag_v(x, y, t, u, v, cbᴰ) = - cbᴰ * v * sqrt(u^2 + v^2)
#bottom_drag_bc_u = FluxBoundaryCondition(bottom_drag_u, field_dependencies=(:u, :v), parameters=cbᴰ)
#bottom_drag_bc_v = FluxBoundaryCondition(bottom_drag_v, field_dependencies=(:u, :v), parameters=cbᴰ)


u_bcs = FieldBoundaryConditions(FluxBoundaryCondition(nothing), top = FluxBoundaryCondition(Q))
v_bcs = FieldBoundaryConditions(FluxBoundaryCondition(nothing))


horizontal_closure = HorizontalScalarDiffusivity(ν=10, κ=10)
# vertical_closure = ScalarDiffusivity(ν=1e-5, κ=1e-5)
vertical_closure = RiBasedVerticalDiffusivity()


model = HydrostaticFreeSurfaceModel(; grid = underlying_grid,
                            coriolis = BetaPlane(latitude = -35),
                            closure = (horizontal_closure,vertical_closure),
                            buoyancy = BuoyancyTracer(),
                            tracers = :b,
                            momentum_advection = WENO(),
                            tracer_advection = WENO(),
                            boundary_conditions = (u=u_bcs, v=v_bcs))

# background density profile based on Argo data
# @inline bg(x,y,z) = (-9.8/1025)*(-0.147 * tanh( 2.6 * ( z + 623 ) / 1000 ) - 1027.6)
# set!(model, b=bg)
set!(model, b=-9.8)


simulation = Simulation(model, Δt = 1minutes, stop_time = 10minutes)
wizard = TimeStepWizard(cfl=1.0, max_change=0.5, max_Δt=3hour)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))


b, = simulation.model.tracers
u, v, w = simulation.model.velocities

extra_outputs = (; 
    u=@at((Center, Center, Center), u),
    v=@at((Center, Center, Center), v),
    w=@at((Center, Center, Center), w),
    Zeta=@at((Center, Center, Center), ∂x(v)-∂y(u)),
)

filename = "../data/output.nc"

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, merge(model.tracers, extra_outputs),
                     overwrite_existing=true,
                     filename = filename,
                     schedule=TimeInterval(1minute))



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(1minute))


run!(simulation)

[simone-silvestri]

As a rule of thumb, you want to use a CFL of 0.2/0.3 with an Adams Bashforth time stepper (the one you are using). A CFL of 1.0 is very large and will lead to numerical instabilities if your system is advection-dominated

[navidcy]

Ok, so that means that Oceananigans does not work for versions 1.8 and earlier, correct? We need to define a compat then. Right @navidcy ?

if we are sure about this then yes

[navidcy]

@iuryt can you strip down the example as much as possible to make a MWE of the problem? I’m finding it a bit hard to distill it. Ideally, a snippet of code that if I paste I will get to see the problem (eg a plot or something).

Question: do you apply wind stress at the boundaries of the domain?

[iuryt]

Of course! I will work on that and post it here. I am only applying the wind stress at the top for zonal velocity.

[navidcy]

Gotcha. Does it vanish at the Bounded boundaries? If not, you expect to see boundary currents there. Although your plots are in longitude so that can’t be the explanation…

[navidcy]

Also, why this value: set!(model, b=-9.8)?

(I'm not saying that it's bad. I just don't see why and I have trouble thinking whether it may be the culprit of blowups...)

[iuryt]

Idk, I think that I was just trying to set a buoyancy for the reference density, thus b = -g rho_0/rho_0 = -g

[navidcy]

Also what you plot are the interpolated velocities at cell centers? If not why are you saving them.

[iuryt]

You are right, I usually do that to help with plotings, but for this case if better to stick to the original face positions.

[navidcy]

OK, from the discussion above, regardless of the boundary layer effect that @iuryt mentions, it also seems that indeed the way we call fkeep! at the moment is broken for Julia v1.8, right? I did change that in #3100 and should have been more careful to run the tests w v1.8 before merging. @simone-silvestri, what's a MWE that would demonstrate this? If that's the case we can either fix it by using something like

if VERSION < v"1.8"
    ...
else
   ...
end

Given that fkeep! is only called on the constructor of the matrix solver. Otherwise we might need to put a compat entry for Julia.

[iuryt]

Hi @navidcy . I will work on a MWE today and post it as an issue, then we can work from that. What do you think?

[navidcy]

You mean for the fkeep? Sounds good.

Did I understood correctly that there is the fkeep issue and also potentially something related with these boundary layer flows you point out? They are two independent things, right?

[iuryt]

I think they are two independent things. The fkeep error occurs when I try to use the bathymetry.

[navidcy]

That’s probably related to the default implicit free surface solver with uneven bathymetry.

[iuryt]

I created #3184 to solve the version problem with fkeep with a minimal reproducible example.

[navidcy]

#3190 deals with the compat entry for Julia v1.9

[iuryt]

Thanks @simone-silvestri! So the plan is to create a function for each velocity at each boundary as suggested by @jagoosw?
Like u[i,1,k],u[i,end,k],v[i,1,k],v[i,end,k],u[1,j,k],u[end,j,k],v[i,1,k], and v[end,j,k]?

[simone-silvestri]

or for the moment you can do

const cbᴰ = 1e-4 # quadratic bottom drag coefficient
@inline lateral_drag_u(x, y, t, u, v, cbᴰ) = - cbᴰ * u * sqrt(u^2 + v^2)
@inline lateral_drag_v(x, y, t, u, v, cbᴰ) = - cbᴰ * v * sqrt(u^2 + v^2)

lateral_drag_bc_u = FluxBoundaryCondition(lateral_drag_u, field_dependencies = (:u, :v), parameters = cbᴰ)
lateral_drag_bc_v = FluxBoundaryCondition(lateral_drag_v, field_dependencies = (:u, :v),  parameters = cbᴰ)

u_bcs = FieldBoundaryConditions(north = lateral_drag_bc_u, south = lateral_drag_bc_u, bottom = FluxBoundaryCondition(nothing), top = FluxBoundaryCondition(Q))
v_bcs = FieldBoundaryConditions(west = lateral_drag_bc_v, east = lateral_drag_bc_v, bottom = FluxBoundaryCondition(nothing), top = FluxBoundaryCondition(nothing))

[glwagner]

The behavior you want (specifying a default bc that is applied to all sides) is not supported with continuous boundary functions, you have to explicitly specify the sides. #3196 should help with the default

Is that a bug? Why isn't that supported?

[iuryt]

Sorry if this is a naive question, but I would imagine that the drag would apply to all velocity components at all applied boundaries. Why north-south only for u and west-east only for v?

[simone-silvestri]

because u is positioned on faces (so exactly on the walls) in the x-direction, so we apply the impenetrability condition (same for v in y-direction)