Replies: 1 comment 5 replies
-
|
To make the But also, Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl Lines 197 to 219 in 1c2a6f8 Computing fluxes across immersed boundaries is definitely a missing feature in the API. Ideally, we would be able to compute this with a 2D kernel, not a 3D kernel. As noted, there will be one flux for |
Beta Was this translation helpful? Give feedback.
5 replies
-
|
Thanks, @glwagner, that was very helpful. I'm still not able to get there either because I'm missing something, or because of a bug somewhere. I'm getting a "method is ambiguous" error. I'll try to work more on it later, but for now, here's a snippet to reproduce the error I'm getting: using Oceananigans
grid_base = RectilinearGrid(topology = (Bounded, Periodic, Bounded),
size = (10, 10, 10),
x = (0, 400), y=(0, 400), z=(0, 100),
halo = (4,4,4))
H = 50; L = 100
mount(x, y) = H * exp(-2((x-grid_base.Lx/2)^2 + (y-grid_base.Ly/2)^2)/L^2)
grid = ImmersedBoundaryGrid(grid_base, GridFittedBottom(mount))
τ = FluxBoundaryCondition(-0.01)
u_bcs = FieldBoundaryConditions(immersed=τ,)
v_bcs = FieldBoundaryConditions(immersed=τ,)
model = NonhydrostaticModel(; grid, closure = SmagorinskyLilly(),
boundary_conditions = (u=u_bcs, v=v_bcs,)
)
set!(model, v=0.05)
u, v, w = model.velocities
fcf = (Face(), Center(), Face())
using Oceananigans.ImmersedBoundaries: conditional_flux, bottom_ib_flux
using Oceananigans.Operators: index_left
using Oceananigans.BoundaryConditions: flip
@inline function conditional_bottom_ib_flux(i, j, k, ibg, bc, loc, c, closure, K, id, clock, fields)
q̃ᴮ = bottom_ib_flux(i, j, k, ibg, bc.bottom, loc, c, closure, K, id, clock, fields)
iᵂ, jˢ, kᴮ = map(index_left, (i, j, k), loc) # Broadcast instead of map causes inference failure
LX, LY, LZ = loc
return conditional_flux(i, j, kᴮ, ibg, LX, LY, flip(LZ), q̃ᴮ, zero(eltype(ibg)))
end
τᵤᶻ_ib = Field(KernelFunctionOperation{Face, Center, Face}(conditional_bottom_ib_flux, grid,
u.boundary_conditions.immersed, fcf, u, model.closure,
model.diffusivity_fields, nothing, model.clock, fields(model)))
compute!(τᵤᶻ_ib)
using Oceananigans.Grids: boundary_node
boundary_node_ccf(i, j, k, grid) = boundary_node(i, j, k, grid, Center(), Center(), Face())
boundary_node_ccf_op = KernelFunctionOperation{Center, Center, Face}(boundary_node_ccf, grid)
bn = Field(boundary_node_ccf_op)
τ = Average(τᵤᶻ_ib, dims=3, condition=boundary_node_ccf_op)And here's the full error: Click me |
Beta Was this translation helpful? Give feedback.
-
|
If we bring the calculation of get_qᵂ(i, j, k, ibg, args...)
get_qᴱ(i, j, k, ibg, args...)
...
@inline function immersed_flux_divergence(i, j, k, ibg::GFIBG, bc, loc, c, closure, K, id, clock, fields)
return div(i, j, k, ibg, loc, get_qᵂ(...), get_qᴱ(...), get_qˢ(...), get_qᴺ(...), get_qᴮ(...), get_qᵀ(...))
end Then it should be relatively easy to implement a user-friendly interface to calculate these fluxes on Oceanostics (although it'd be a 3D kernel for now). The thing is that I think the If you agree with this let me know and I'll open a PR. |
Beta Was this translation helpful? Give feedback.
-
|
I think that method ambiguity is a bug. @simone-silvestri might know more. I wonder if that is tested. |
Beta Was this translation helpful? Give feedback.
-
|
Hmmm indeed it looks like a bug, I ll take a look tomorrow. |
Beta Was this translation helpful? Give feedback.
-
|
The proper implementation of this should launch a 2D kernel that knows the |
Beta Was this translation helpful? Give feedback.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
I have a simulation that can be illustrated with the short code below, where I'm setting up a Gaussian bump in the middle and starting with a flow of
0.01m/sin theydirection:I'm trying to output the drag force (really the drag work is what I want, but let's stick with the force for now) induced by this object. Starting with the z-direction u-drag, I thought I could do something like:
which (I think) should isolate the value of the fluxes only at the boundary of my immersed solid.
However, the code above doesn't work and I think it's because
conditionactually has to be anAbstractOperation(is that right?). (Error message available below:) If so, then I'm not sure how to do the above.Click to expand
Furthermore, since the topography above is a symmetric Gaussian bump, it's got x and y direction fluxes on both sides of x and y, which would require two separate outputs. I'm really not sure how to deal with that tbh.
So, can I make this method work? And maybe a better question: is there a better way to accomplish what I want?
Beta Was this translation helpful? Give feedback.
All reactions