-
Notifications
You must be signed in to change notification settings - Fork 204
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Adds a pertubation advection open boundary matching scheme #3977
base: main
Are you sure you want to change the base?
Conversation
@tomchor, do you think something like this: Oceananigans.jl/validation/open_boundaries/cylinder.jl Lines 18 to 76 in 32cd354
could belong in Oceanostics? |
I have also run the same case with different domain lengths (6m, 12m - shown, 18m and 30m) and they all produce very similar Cd and St |
Sure! Why not? :) |
@inline function _fill_east_halo!(j, k, grid, u, bc::PAOBC, loc::Tuple{Face, Any, Any}, clock, model_fields) | ||
i = grid.Nx + 1 | ||
|
||
Δt = clock.last_stage_Δt | ||
|
||
Δt = ifelse(isinf(Δt), 0, Δt) | ||
|
||
Δx = xspacing(i, j, k, grid, loc...) | ||
|
||
ūⁿ⁺¹ = getbc(bc, j, k, grid, clock, model_fields) | ||
|
||
uᵢⁿ = @inbounds u[i, j, k] | ||
uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k] | ||
|
||
U = max(0, min(1, Δt / Δx * ūⁿ⁺¹)) | ||
|
||
τ = ifelse(ūⁿ⁺¹ >= 0, | ||
bc.classification.matching_scheme.outflow_timescale, | ||
bc.classification.matching_scheme.inflow_timescale) | ||
|
||
|
||
τ̃ = Δt / τ | ||
|
||
uᵢⁿ⁺¹ = (uᵢⁿ + U * uᵢ₋₁ⁿ⁺¹ + ūⁿ⁺¹ * τ̃) / (1 + τ̃ + U) | ||
|
||
@inbounds u[i, j, k] = uᵢⁿ⁺¹#ifelse(active_cell(i, j, k, grid), uᵢⁿ⁺¹, zero(grid)) | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we could implement everything using getindex()
so that we implement it once (or at max twice; one for left and other for right BC) and can use it in three dimensions.
For example this
uᵢⁿ = @inbounds u[i, j, k]
uᵢ₋₁ⁿ⁺¹ = @inbounds u[i - 1, j, k]
would become
uᵢⁿ = @inbounds getindex(u, boundary_indices...)
uᵢ₋₁ⁿ⁺¹ = @inbounds getindex(u, one_off_boundary_indices...)
This would avoid errors when transcribing to other dimensions (I remember catching a couple for the flat extrapolation matching scheme).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yeah, I was trying to think of a clean way to do this
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Attempted this now (but I've still not written out the other directions in case we want to change things first)
0e645c5
to
32cd354
Compare
afecfb3
to
e66f3d1
Compare
There's a typo in your top equation right? I think it should be |
What do you mean? What is the consequence of this approximation? |
I was thinking about this more and really to remain valid then in : as I found a typo in the west boundary which might be why the oscillating case wasn't working, rerunning now |
|
Okay, I agree with you. You can also analyze the update formula / backward Euler step in the limit U -> 0. It doesn't look like it needs to be regularized in any way in that case though, it looks ok. |
Is there a reference for this method? If there is, please include it on the methods docstring. If not, it'd be good to include a brief explanation like the one at the top comment in the dosctring. I don't think this radiation method is very well known. |
|
||
u = normal_velocity(Val(orientation), model) | ||
|
||
CUDA.@allowscalar u.data .= 1 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we're trying to avoid introducing new @alloscalar
instances to testing. @glwagner can you confirm?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought the problem was that tests currently have a global @allowscalar
, but some tests need them or will be more complicated to write?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you can write most of them simply without allowscalar
, but I'm not an expert on this. For example you may be able to write something like interior(u) .= 1
without it (I haven't tested it though), and it's a pretty simple change.
I say try to avoid it with simple statements like the one I just mentioned. If we can't do it simply and @glwagner doesn't manifest any opinions on it, we move forward. How does that sound?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I sorted this now, it wasn't as straight forward as just using interior because I needed to set a peripheral node, but I just used view
instead
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Made a few comments, but it looks like there are old unresolved comments and some files here and there which do not seem to be ready for prime yet. let me know if anything needs clarification.
Co-authored-by: Tomas Chor <[email protected]>
@tomchor I think I've addressed your comments now |
|
||
ax = Axis(fig[1, 1], aspect = DataAspect()) | ||
simulation.output_writers[:drag] = JLD2OutputWriter(model, (; drag_force), | ||
schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
schedule = IterationInterval(Int(0.1/Δt)),#TimeInterval(0.1), | |
schedule = IterationInterval(Int(0.1/Δt)), |
outputs = (; u, v, p, ζ) | ||
|
||
simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs, | ||
schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
schedule = IterationInterval(Int(2/Δt)),#TimeInterval(0.1), | |
schedule = IterationInterval(Int(2/Δt)), |
|
||
progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))" | ||
simulation = Simulation(model; Δt, stop_time) | ||
#conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
#conjure_time_step_wizard!(simulation, cfl=1.0, IterationInterval(3); max_Δt) |
# there is some problem using ConjugateGradientPoissonSolver with TimeInterval because the timestep can go really small | ||
# so while I identify the issue I'm using IterationInterval and a fixed timestep |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Have you tried setting minimum_relative_step
? (see #3606)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This seems to fix the issue! I hadn't tried this since early november/december so I hadn't seen you'd made this fix so this is great, thank you
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll update the validation script
Co-authored-by: Tomas Chor <[email protected]>
We begin with the equation governing the fluid in the interior: | ||
∂ₜu⃗ + u⃗⋅∇u⃗ = −∇P + F⃗, | ||
and note that on the boundary the pressure gradient is zero. | ||
We can then assume that the flow composes of mean (U⃗) and pertubation (u⃗′) components, | ||
and considering the x-component of velocity, we can rewrite the equation as | ||
∂ₜu₁ = -u₁∂₁u - u₂∂₂u₁ - u₃∂₃u₁ + F₁ ≈ - U₁∂₁u₁′ - U₂∂₂u₁′ - U₃∂₃u₁′ + F. | ||
|
||
Simplify by assuming that U⃗ = Ux̂, an then take a numerical step to find u₁. | ||
|
||
When the boundaries are filled the interior is at time tₙ₊₁ so we can take | ||
a backwards euler step (in the case that the mean flow is boundary normal) on a right boundary: | ||
(Uⁿ⁺¹ - Uⁿ) / Δt + (u′ⁿ⁺¹ - u′ⁿ) / Δt = - Uⁿ⁺¹ (u′ⁿ⁺¹ᵢ - u′ⁿ⁺¹ᵢ₋₁) / Δx + Fᵤ. | ||
|
||
This can not be solved for general forcing, but if we assume the dominant forcing is | ||
relaxation to the mean velocity (i.e. u′→0) then Fᵤ = -u′ / τ then we can find u′ⁿ⁺¹: | ||
u′ⁿ⁺¹ = (uⁿ + Ũu′ⁿ⁺¹ᵢ₋₁ - Uⁿ⁺¹) / (1 + Ũ + Δt/τ), | ||
|
||
where Ũ = U Δt / Δx, then uⁿ⁺¹ is: | ||
uⁿ⁺¹ = (uᵢⁿ + Ũuᵢ₋₁ⁿ⁺¹ + Uⁿ⁺¹τ̃) / (1 + τ̃ + U) | ||
|
||
where τ̃ = Δt/τ. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this not rendered correctly too @tomchor ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If I use FiraCode on chrome then some of the stuff in the first line of code doesn't render correctly:
On firefox it all renders fine:
And when I manually switch to JuliaMono it all renders fine.
But after finding out that this is all browser/machine dependent, I do think we shouldn't spend too much time trying to "fix" this. I say we move on from tweaking the unicode and address the rest of the comments so that we can merge.
This PR adds a
PertubationAdvection
open boundary condition which assumes that there is some mean flow (i.e. prescribed in a channel or from a coarser resolution model) which is homogeneous on the boundary but can vary in time.On the boundary, we approximate (for an x-boundary) to:
where$u=u\prime+U$
I have chosen to take an explicit/backwards Euler step:
because then we don't have to store any other information (i.e. if we did a forward step we would need$u\prime(x_{b-1}, t^n)$ ). This results in:
where$\tilde{U}=\min\left(1, U\Delta t / \Delta x\right)$ . I have also added restoring to $U$ (i.e. damping $u\prime\to0$ ) in some timescale $\tau$ which can be different for inflow and outflow to allow the scheme to work when there is inflow, which results in:
where$\tilde{\tau}=\Delta t / \tau$ .
This has the limitation that if we want to extend to have wall tangent mean velocity ($\partial_t u \approx -(\vec{U}\cdot\nabla) u$ ) then we either have to solve the whole boundary at once since every point is going to depend on its boundary neighbours, or we need to take a forward Euler step in which case we need the first interior point at the previous time step.
I think we probably will need to address this for real nesting cases at some point.
u.mp4
For now, it seems to be working okay (left plot vorticity, right plot x-velocity):
(old plot)
https://github.com/user-attachments/assets/21ae3eb2-5d3b-4c08-86fb-bb9c222e8e6aI have also checked the Strouhal number which matches exactly to the expected values (~0.17).