Velocity profile of channel flow case has difference with log-law

[Tinydog8]

Hi all, I have met some strange things in a simple channel flow case, the velocity profile is larger than the log-low profile, and the momentum flux of the second and third points are obviously lower than bottom boundary condition (or other point near the bottom boundary). The code is written below, what causes this difference?

const H=15 #/m
grid = RectilinearGrid(GPU(),size=(64,64,64), extent=(π*H, π*H, H))
const u★=0.01 #friction velocity
Fx(x,y,z,t)=u★^2/H #forcing

const z₀ = H*1e-4 # m (roughness length)
const κ = 0.4 # von Karman constant
const z₁ = -1*znodes(Center,grid)[grid.Nz] # Closest grid center to the bottom
const cᴰᵇ = (κ / log(z₁ / z₀))^2 # Drag coefficient

@inline drag_u(x, y, t, u, v, p) = - p.cᴰᵇ * √(u^2 + v^2) * (u)
@inline drag_v(x, y, t, u, v, p) = - p.cᴰᵇ * √(u^2 + v^2) * (v)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰᵇ))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰᵇ))

u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(0.0),bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

model = NonhydrostaticModel(; grid, coriolis,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers =(:T,:S),
                            buoyancy = SeawaterBuoyancy(),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs,v=v_bcs,T=T_bcs,S=S_bcs),
                            forcing=(u=Fx,))

[simone-silvestri]

Hello, Do you want to simulate a channel flow or a turbulent boundary layer? Because at the moment, you are just specifying the drag at the bottom (which I guess is a model for the no-slip boundary condition), while a channel requires no-slip at both the top and the bottom. Also, in case you want to simulate a turbulent boundary layer, take care that it is not a periodic case but spatially developing in the streamwise direction, so if you are looking for the self-similar solution (the log-law), you need to rescale somewhere at the beginning or the end of the domain. This is not a problem with a fully developed channel flow that is statistically homogeneous in the streamwise direction.

Given that you probably need to specify the drag at the top and bottom (unless drag is only for roughness), also cᴰᵇ is wrong. This is because the closest grid center to the bottom is the first element in the znodes array, not the last:

const z₁ = -1*znodes(Center,grid)[1] # Closest grid center to the bottom

Let me know if that works

[Tinydog8]

Thank you for your reply! And apologies for have not been clear. Firstly, we want to simulate a "half-channel flow", so we only use no-slip boundary condition at the bottom. Secondly, we have used periodic boundary condition in x (streamwise) and y (crosswise) direction. Finally, in this case the z-axis is range from -15 to 0, so I can't use const z₁ = -1*znodes(Center,grid)[1] because if I do that, the value of z1 will change to 14.88. My code const z₁ = -1*znodes(Center,grid)[grid.Nz] will keep z1 equals to Closest grid center to the bottom (0.12). Everything seems right, but the result seems strange!

[chabbymark]

Hi all, I have met some strange things in a simple channel flow case, the velocity profile is larger than the log-low profile, and the momentum flux of the second and third points are obviously lower than bottom boundary condition (or other point near the bottom boundary). The code is written below, what causes this difference?

const H=15 #/m
grid = RectilinearGrid(GPU(),size=(64,64,64), extent=(π*H, π*H, H))
const u★=0.01 #friction velocity
Fx(x,y,z,t)=u★^2/H #forcing

const z₀ = H*1e-4 # m (roughness length)
const κ = 0.4 # von Karman constant
const z₁ = -1*znodes(Center,grid)[grid.Nz] # Closest grid center to the bottom
const cᴰᵇ = (κ / log(z₁ / z₀))^2 # Drag coefficient

@inline drag_u(x, y, t, u, v, p) = - p.cᴰᵇ * √(u^2 + v^2) * (u)
@inline drag_v(x, y, t, u, v, p) = - p.cᴰᵇ * √(u^2 + v^2) * (v)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰᵇ))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰᵇ))

u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(0.0),bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

model = NonhydrostaticModel(; grid, coriolis,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers =(:T,:S),
                            buoyancy = SeawaterBuoyancy(),
                            closure = AnisotropicMinimumDissipation(),
                            boundary_conditions = (u=u_bcs,v=v_bcs,T=T_bcs,S=S_bcs),
                            forcing=(u=Fx,))

I also have the same problem here. I am trying to simulate the neutral turbulent boundary layer here. However, I found that the velocity shear at the first grid points is much larger than that predicted by the Monin-Obukhov similarity theory.

Any ideas? @glwagner @tomchor

[tomchor]

Here are a couple of comments that come to mind when looking at your code.

1. You're using a call to znodes(Center,grid), which indicates that you're using an out-of-date version of the code. I suggest you update the code to the latest version and then try again. It's good to always keep your code up-to-date (especially when posting here) because the code is always being improved upon (and in some instances some bugs get fixed which may be important to your example!)

2. That said, about this:

Finally, in this case the z-axis is range from -15 to 0, so I can't use const z₁ = -1*znodes(Center,grid)[1] because if I do that, the value of z1 will change to 14.88. My code const z₁ = -1*znodes(Center,grid)[grid.Nz] will keep z1 equals to Closest grid center to the bottom (0.12).

This isn't the best way to define $z_1$. The current statement will break if you, for example, decide to stretch the grid vertically. The more robust way to define $z_1$ is

julia> grid = RectilinearGrid(CPU(), size=(64,64,64), extent=(π*H, π*H, H))
64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [1.54466e-16, 47.1239) regularly spaced with Δx=0.736311
├── Periodic y ∈ [1.54466e-16, 47.1239) regularly spaced with Δy=0.736311
└── Bounded  z ∈ [-15.0, 0.0]           regularly spaced with Δz=0.234375

julia> using Oceananigans.Grids: zspacing

julia> zspacing(1, 1, 1, grid, Center(), Center(), Center())/2 # Half the distance around a center point
0.1171875

This is clearer, and will give you the correct value even if you change the grid in the future. Note that in the call to zspacing() above, we need to specify the x, y and z positions to get the spacing, but in a RectilinearGrid Δz is constant in x and y so those arguments don't really matter.

Also note that the code above is written for the most recent Oceananigans version (for which you'll need Julia 1.9).

3. Now to your main point (sorry for the tangents! :grimacing:), I don't see anything obviously wrong with the code. Just from looking at it, it looks like it should work. (Btw, when I say that, I assume this is a only snippet, and not the full-code, since I get errors when I try to run your code (coriolis, T_bcs and S_bcs aren't defined, Fx looks like it'd cause an error, etc.), so I haven't been able to properly run and test your snippet.)

That said, there are things I personally would try:

• Have you tried better resolving the model? A 64^3 model should be able to achieve a log-law, but it's still worth investigating.
• Have you tried using the SmagorinskyLilly closure? I've had some issues with AMD in the past (although not exactly the same as yours), and I remember @glwagner mentioning having an issue with AMD near boundaries that could explain the behavior you're seeing.
• Finally, I'd suggest starting with a very simple model first for debugging/investigating. Just the bare minimum to reproduce a log-law. (i.e., no tracers, buoyancy, Coriolis, forcing...)

@Tinydog8 @chabbymark Hope this helps!

[glwagner]

However, I found that the velocity shear at the first grid points is much larger than that predicted by the Monin-Obukhov similarity theory.

I don't have much experience with solid wall boundary layers.

A quick search returns this review:

https://journals.aps.org/prfluids/pdf/10.1103/PhysRevFluids.2.104601?casa_token=VBO0hrNqp-cAAAAA%3A3XEKZfLTdhiVluwRl8sCQCkOK44GoGUX-546uFtqQjSWAUIDQpKyyQsA4lQd65Oz6Kw5ClWias_CVQ0

suggesting that correct flux predictions in wall-modeled LES is unsolved. For example:

shows a mismatch between DNS and wall-modeled LES for a few standard codes.

Based on the literature, a failure to reproduce the log-law is expected?

It might help if you provide more background on what exactly you're trying to achieve, and why you believe the current approach will achieve that objective. For example, are we attempting to reproduce a known result?

[chabbymark]

@glwagner @tomchor I have just tested the turbulent boundary layer with wall model using Oceananigans on Julia 1.6.7 and Julia 1.9.2. The performances of AMD in both versions are exactly the same. It turns to overestimate the velocity shear at the second node from the wall, therefore, turns to overestimate the velocity in the middle and upper parts of the boundary layers.

I used to use the SGS model based on the Lagrangian-averaged scale-dependent dynamic model (LASD) (https://pubs.aip.org/aip/pof/article/17/2/025105/895722/A-scale-dependent-Lagrangian-dynamic-model-for). @tomchor is very familiar with this SGS model. The performance of the LASD close to the wall is usually good, as you can see here .

I guess the problem of AMD is partly solved in reference Yang et al. (2017). Now the problem is that if someone can implement this filtering in the code or not. I am stilling learning the Oceananigans and Julia. I hope that someday in the future, I am able to implement this technique.

[glwagner]

Hmm ok. But just to make sure I understand --- you are saying that you disagree with Yang, Park, and Moin (2017) paper, correct? In that paper, they apparently use the Lagrangian-scale-averaged (LSA) SGS model

Of course, whether or not this LSA closure has an issue, there could still be a problem with wall-modeled stresses when using AnisotropicMinimumDissipation. (In the above paper, they claim that the SGS model doesn't matter.)

I'm curious what happens if you try

advection = WENO(order=9)

You may also need to add something like halo = (4, 4, 4) or halo = (5, 5, 5) to the RectilinearGrid constructor.

I'm also curious what happens if you use advection = WENO(order=9) with closure = nothing.

Now the problem is that if someone can implement this filtering in the code or not.

Are you referring to the time-filtering described in equation 2?

It is not necessary to change Oceananigans to implement this --- it can be implemented in your script. You need to introduce a new field to store the time-filtered wall velocity, and update this field in a Simulation.callbacks. It should be straightforward.

By the way, I think this should be a Discussion (rather than an issue), because this does not (as far as I can tell) involve a bug of some kind in Oceananigans. It could motivate implementing a new SGS model, but that's a discussion for another day perhaps...

[glwagner]

@Tinydog8 let me know if you're ok to convert this to discussion.

Here's a sketch of how to accumulate u_wm:

u_wall_model = XFaceField(grid, indices=(:, :, 1)) # It's purely a niceity to specify `indices` properly here --- the only thing we really need is a 2D array

# Also define other velocity components as needed
# Define drag boundary condition in terms of u_wall_model, other velocity components, etc using the `discrete_form=true`

# This function is meant to be used in a callback
function compute_wall_model_velocity!(sim)
    u, v, w = sim.model.velocities

    filtering_time_scale = 1.0 # seconds --- for example
    ϵ = sim.Δt / filtering_time_scale

    u_wm = interior(u_wall_model, :, :, 1) # extract a view to broadcast over
    u_LES = interior(u, :, :, 1) # a view into the LES velocity within the first grid cell
    @. u_wm = (1 - ϵ) * u_wm + ϵ * u_LES

    # Compute other components as needed
    return nothing
end

simulation.callbacks[:wall_model] = Callback(compute_wall_model_velocity!)

An even more slick approach would build new auxiliary fields for the wall model velocities. In that case (if I'm not mistaken), then the wall model velocities are accessible from the argument to the boundary condition function (and don't have to be referenced as global variables).

[tomchor]

@glwagner if this works (i.e. if AMD with this change indeed reproduces the log-law better) this would make a pretty cool example for the docs.

[glwagner]

@glwagner if this works (i.e. if AMD with this change indeed reproduces the log-law better) this would make a pretty cool example for the docs.

We can put it on the list though I'd be hesitant to move too quickly because our docs will probably be getting a lot heavier in the near future with sphere examples, and examples with more complex bathymetry. Perhaps you can make the tilted bottom boundary layer 3D and add it into that one, so we don't pay the price of a new independent example?

[glwagner]

Is this still an issue? Should we convert this to a discussion?

[Tinydog8]

Is this still an issue? Should we convert this to a discussion?

Sorry for the delay! Yes, that is still a problem and you can convert this to a discussion. I will try to use higher order advcetion scheme later. Thanks a lot!