Bottom Drag Boundary Conditions

[IsaiahCuadras]

Hi Everyone!

I am just getting started with Oceananigans and I am trying to run an LES case of a stratified bottom Ekman boundary layer. I am doing this by making some modifications to the Langmuir turbulence example and the tilted bottom boundary layer example.

I believe that I am having some trouble with implementing the bottom drag in particular. I actually can't tell if the issue is with my implementation of the model, or if there is just an issue in the way I am loading the data. The code below is what I am using to build the simulation (omitting initial conditions and outputting functions from this post for now, but can include if needed):

using Oceananigans
using Oceananigans.Units: minute, minutes, hours
using Printf

# Parameters 
coriolis = FPlane(f=1e-4) # s⁻¹
N² = 1.936e-5 
u_bckgnd = 0.06
v_max = 0.01

# Set grid
grid = RectilinearGrid(size=(32, 32, 32), extent=(128, 128, 64), topology=(Periodic, Periodic, Bounded))

# Boundary Conditions
z₀ = 0.1 # m (roughness length)
κ = 0.4  # von Karman constant
z₁ = abs(first(znodes(grid, Center()))) # Closest grid center to the bottom
cᴰ = (κ / log(z₁ / z₀))^2 # Drag coefficient
@inline drag_u(x, y, t, u, v, p) = - p.c_drag * √(u^2 + v^2) * u
@inline drag_v(x, y, t, u, v, p) = - p.c_drag * √(u^2 + v^2) * v
drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(c_drag = cᴰ, u_inf = u_bckgnd))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(c_drag = cᴰ, u_inf = u_bckgnd))

vel_bc_u = FieldBoundaryConditions(bottom=drag_bc_u, top=GradientBoundaryCondition(0.0))
vel_bc_v = FieldBoundaryConditions(bottom=drag_bc_v, top=GradientBoundaryCondition(0.0))
buoy_bcs = FieldBoundaryConditions(top = GradientBoundaryCondition(0),
                                bottom = FluxBoundaryCondition(0))
                                
# Background flow forcing
@inline v_forcing_func(x, y, z, t, p) = p.u_inf*p.cor
v_forcing = Forcing(v_forcing_func, parameters=(cor = coriolis.f, u_inf = u_bckgnd))

# Instantiate Model
model = NonhydrostaticModel(; grid, coriolis,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=vel_bc_u, v=vel_bc_v, b=buoy_bcs),
                            forcing = (v=v_forcing,))

The model runs without error, but trying to follow the Langmuir turbulence example for post-processing returns the following error at the following line of code:

vₙ = @lift time_series.v[$n]

ERROR: MethodError: no method matching FieldBoundaryConditions(::Tuple{Colon, UnitRange{Int64}, Colon}, ::Missing)

However, if I change the boundary condition to a slip or no-slip boundary condition, this error goes away and the post-processing script runs fine. Any idea why this might be?

Apologies for what I am sure is a silly question!

Best,

Isaiah

[glwagner]

Not silly! I think this is a bug. This might fix it: #4557

[navidcy]

From what I understand there is another error that @IsaiahCuadras is pointing to which comes from view. I get:

julia> view(time_series.u[1], :, :, 10)
ERROR: MethodError: no method matching FieldBoundaryConditions(::Tuple{Colon, Colon, UnitRange{Int64}}, ::Missing)

Closest candidates are:
  FieldBoundaryConditions(::Tuple, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any)
   @ Oceananigans ~/Library/CloudStorage/OneDrive-TheUniversityofMelbourne/Documents/Research/Oceananigans.jl-v3/src/BoundaryConditions/field_boundary_conditions.jl:60
  FieldBoundaryConditions(::W, ::E, ::S, ::N, ::B, ::T, ::I) where {W, E, S, N, B, T, I}
   @ Oceananigans ~/Library/CloudStorage/OneDrive-TheUniversityofMelbourne/Documents/Research/Oceananigans.jl-v3/src/BoundaryConditions/field_boundary_conditions.jl:51
  FieldBoundaryConditions(::MultiRegionGrid, ::Any, ::Any; kwargs...)
   @ Oceananigans ~/Library/CloudStorage/OneDrive-TheUniversityofMelbourne/Documents/Research/Oceananigans.jl-v3/src/MultiRegion/multi_region_field.jl:154
  ...

Stacktrace:
 [1] view(f::Field{…}, i::Function, j::Function, k::Int64)
   @ Oceananigans.Fields ~/Library/CloudStorage/OneDrive-TheUniversityofMelbourne/Documents/Research/Oceananigans.jl-v3/src/Fields/field.jl:346
 [2] top-level scope
   @ REPL[86]:1
Some type information was truncated. Use `show(err)` to see complete types.

But after #4547 I get:

julia> view(time_series.u[1], :, :, 10)
32×32×1 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 32×32×32 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: Nothing 
├── indices: (:, :, 10:10)
└── data: 38×38×1 OffsetArray(view(::Array{Float64, 4}, :, :, 13:13, 1), -2:35, -2:35, 10:10) with eltype Float64 with indices -2:35×-2:35×10:10
    └── max=0.0595074, min=0.0594558, mean=0.059481

[IsaiahCuadras]

Just want to confirm that @navidcy's most recent remark is indeed the case. I had not come across this since my script was stopped by the original error mentioned in this post (now fixed by PR #4545).

@navidcy is your second piece of code indicating that the view error has been fixed as well? I updated to the most recent oceananigans version, but I still get that MethodError from view.

[navidcy]

It’s fixed on main branch. I can tag a new release tomorrow

[navidcy]

https://github.com/CliMA/Oceananigans.jl/releases/tag/v0.96.30

Should have fixed it!

[IsaiahCuadras]

Yep, that did the trick. Thanks so much!

[glwagner]

Not sure if #4545 is related, @navidcy