Resolved LES in Oceananigans

[mmr0]

I'm after some help setting up a resolved LES in Oceananigans. By resolved LES i'm specifically referring to resolving the viscous sublayer/wall stress for a wall-bounded flow, e.g. no wall model needed. Has anyone done this, and if so are you willing to share your experience?

Where I'm at: I'm currently trying to replicate a DIABLO simulation of an upside down bottom Ekman layer (i.e. the no-slip boundary is at the top of the domain) as a test but despite having the same resolution as I've used with DIABLO I don't have anything that looks like a viscous sublayer in my simulation, and the stress at the interface is also not what I expect. This was using closure = AnisotropicMinimumDissipation(). Code is below if interested. Note that the code runs and i get pretty sensible looking outer boundary layer and mixed layer, it's just the near-wall stuff that I'm a bit dubious of.

A potentially related issue: I was looking into the Tuple closure method, which would suit me since I'd like to modify the molecular viscosity to something suitable for Polar oceans. I tried:
amd_closure = (AnisotropicMinimumDissipation(), ScalarDiffusivity(ν=2.0e-6, κ=1.4e-7))
but got the error type Tuple has no field νₑ. I saw that this was likely associated with a bug fixed a few versions ago by @tomcho, so I updated updated to the latest version, but I still get the error.

A further potentially related issue: Also came across the issue #1879 regarding tests that the tuple closure method works, and wondered if these were ever carried out?

Sorry for the lengthy post! This is my first foray into Oceananigans, so I'm really not sure what's just user error at this point. Cheers!

using Random
using Printf
using CairoMakie
using Oceananigans
using Oceananigans.Units: minute, minutes, hour

Nz = 300          # number of points in the vertical direction
Lz = 6          # (m) domain depth

# Stretched grid using tanh
h(k) = (k - 1) / Nz
const CSY = 2.5
z_faces(k) = Lz * (tanh(CSY * h(k)) / tanh(CSY) - 1)

grid = RectilinearGrid(GPU();
                       size = (128, 128, Nz),
                          x = (0, 6),
                          y = (0, 6),
                          z = z_faces)

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = 2e-4,
                                                                    haline_contraction = 8e-4))
dTdz = 0.0015 # K m⁻¹
T_bcs = FieldBoundaryConditions(top = GradientBoundaryCondition(0.0),
                                bottom = GradientBoundaryCondition(dTdz))
u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(0.0))
v_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(0.0))
w_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(0.0))
S_bcs = FieldBoundaryConditions(top=GradientBoundaryCondition(0.0))

const u_m = 0.014
coriolis = FPlane(f=1.4e-4)

@inline geostrophy(x, y, z, t, p) = @inbounds -p.U0 * p.f
v_forcing = Forcing(geostrophy, parameters = (U0=-u_m, f=coriolis.f))

amd_closure = (AnisotropicMinimumDissipation(), ScalarDiffusivity(ν=2.0e-6, κ=1.4e-7))

model = NonhydrostaticModel(; grid, buoyancy,
                            advection = advection = UpwindBiasedFifthOrder(),
			    timestepper = :RungeKutta3,
                            tracers = (:T, :S),
                            coriolis = coriolis,
                            closure = amd_closure, #AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs,v=v_bcs, T=T_bcs, S=S_bcs),
			    forcing = (v=v_forcing,)
			    )

# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

# Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = 20 + dTdz * z + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

# Velocity initial condition - needs to satisfy u_bc and gesostrophy
u_init(z) = -u_m * tanh(z*3)

# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = u_init(z) + u_m*1e-3 * Ξ(z)
wᵢ(x, y, z) = u_m * 1e-3 * Ξ(z)
# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=35)

[glwagner]

We believe this use case is important and so we explicitly test it:

Oceananigans.jl/test/test_turbulence_closures.jl

Lines 137 to 145 in 2fb0f77

	function time_step_with_tupled_closure(FT, arch)
	closure_tuple = (AnisotropicMinimumDissipation(FT), ScalarDiffusivity(FT))
	model = NonhydrostaticModel(closure=closure_tuple,
	grid=RectilinearGrid(arch, FT, size=(1, 1, 1), extent=(1, 2, 3)))
	time_step!(model, 1, euler=true)
	return true
	end

in addition to some other tests for tuple closures. That said, there is probably room to expand the tests for tuple closures. But I'm not sure the tuple itself is your issue.

I tried: amd_closure = (AnisotropicMinimumDissipation(), ScalarDiffusivity(ν=2.0e-6, κ=1.4e-7)) but got the error type Tuple has no field νₑ.

This error is probably associated with the code you are using for diagnostics. Can you share that?

[glwagner]

Nice. Since the numerics are different, maybe we also need to test multiple resolutions?

I'd be curious about the effect of different combinations of LES closures and advection schemes. For example, we usually run LES with advection = WENO() and the AnisotropicMinimumDissipation() closure. However, it would be also interesting to test advection = CenteredSecondOrder() (arguably closer to what DIABLO uses), and also advection = WENO() with no LES closure (ie just closure = ScalarDiffusivity()).

[glwagner]

Are you running with UpwindBiasedFifthOrder()? If so I recommend switching to WENO as a first step.

[mmr0]

Yes I recently switched to WENO on @zhazorken recommendation! Looks similar at first glance. I will investigate some different combinations and report back!

[glwagner]

Another thing that might be important if you have very fine grid spacing near the wall is VerticallyImplicitTimeDiscretization:

vitd = VerticallyImplicitTimeDiscretization()
molecular_diffusivity = ScalarDiffusivity(vitd, ...)

you may also want to use for AMD.

[mmr0]

Thanks @glwagner - I am just coming back to this now so great timing!

[navidcy]

Oh I just saw this! Nice to see you started running Oceananigans @mmr0 !