Skip to content
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.

Sign up for GitHub

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

Allow divergent Lagrangian-mean velocity by subtracting Stokes drift before pressure correction #4013

Draft
wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Allow divergent Lagrangian-mean velocity by subtracting Stokes drift before pressure correction #4013

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0fc5708
Allow divergent Lagrangian-mean velocity by subtracting Stokes drift …
glwagner Dec 18, 2024
b95f885
Update StokesDrifts.jl
BrodiePearson Dec 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 39 additions & 3 deletions src/Models/NonhydrostaticModels/pressure_correction.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,21 +1,57 @@
import Oceananigans.TimeSteppers: calculate_pressure_correction!, pressure_correct_velocities!

const c = Center()
const f = Face()

using Oceananigans.Grids: node
using Oceananigans.StokesDrifts: StokesDrift, parameters_tuple

function modulate_by_stokes_drift!(model, sgn, t=time(model))
stokes_drift = model.stokes_drift
if stokes_drift isa StokesDrift
grid = model.grid
arch = architecture(grid)
u, v, w = model.velocities
launch!(arch, grid, :xyz, _modulate_by_stokes_drift, u, v, w, sgn, grid, t, stokes_drift)
end
return nothing
end

@kernel function _modulate_by_stokes_drift(u, v, w, sgn, grid, time, sd)
i, j, k = @index(Global, NTuple)

pt = parameters_tuple(sd)
Xu = node(i, j, k, grid, f, c, c)
Xv = node(i, j, k, grid, c, f, c)
Xw = node(i, j, k, grid, c, c, f)

@inbounds begin
u[i, j, k] = u[i, j, k] + sgn * sd.uˢ(Xu..., time, pt...)
v[i, j, k] = v[i, j, k] + sgn * sd.vˢ(Xv..., time, pt...)
w[i, j, k] = w[i, j, k] + sgn * sd.wˢ(Xw..., time, pt...)
Comment on lines +29 to +31
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this subtracting/adding the Stokes drift at a specific depth from a depth-averaged current? Perhaps this is the kind of issue you were alluding to @qingli411

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think here we should use cell-averaged Stokes drift and in the finite volume approximation of $\partial_z u^S$ we will need the values at specific locations. When Stokes drift is specified as an analytical function, we can make sure this is done consistently. But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

Copy link
Member Author

@glwagner glwagner Dec 22, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, it's a good idea to discuss the discretization. I think on a C-grid we probably want to have the Stokes drift at the same location as the velocity components. This would not be cell centers, but Face, Center, Center for u^S, Center, Face, Center for v^S, etc.

Probably, there is a correct discretization of the term $\tilde \omega \times u^L$ that conserves energy (where $\tilde \omega = \nabla \times u^S$ is the curl of the Stokes drift), similar to how there is a correct discretization of vorticity / Coriolis that conserves energy.

For example see how we do it for Coriolis on the sphere:

Oceananigans.jl/src/Coriolis/hydrostatic_spherical_coriolis.jl

Lines 112 to 119 in d8273e7

@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
@inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)
@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
@inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)

If we want to be able to naturally compute $\omega$ at the "C-grid vorticity locations" then I believe we want u^S at the same location as the velocity field? The three components of vorticity are at cff, fcf, and ffc for the x, y, z components of vorticity respectively.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can build a simple example to demonstrate conservation of energy (which is $(u^L)^2 / 2$) which will show that we've done the discretization correctly.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

That's a fair point. Though it would still be possible to compute the Stokes drift numerically by integrating its derivatives, but maybe that's not so convenient.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

That's a fair point. Though it would still be possible to compute the Stokes drift numerically by integrating its derivatives, but maybe that's not so convenient.

The pressure solver depends only on the divergence of the velocity field. If the Stokes drift gradients were provided, would we now need to specify every derivative since the terms needed for the pressure solver Eulerian-ization ($\partial_xu^s$, $\partial_yv^s$ and $\partial_zw^s$) are distinct from the ones needed for the wave-averaged forces:

Oceananigans.jl/src/StokesDrifts.jl

Lines 153 to 167 in 0fc5708

struct StokesDrift{P, US, VS, WS, VX, WX, UY, WY, UZ, VZ, UT, VT, WT}
uˢ :: US
vˢ :: VS
wˢ :: WS
∂x_vˢ :: VX
∂x_wˢ :: WX
∂y_uˢ :: UY
∂y_wˢ :: WY
∂z_uˢ :: UZ
∂z_vˢ :: VZ
∂t_uˢ :: UT
∂t_vˢ :: VT
∂t_wˢ :: WT
parameters :: P
end

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK. Maybe asking the users to specify every derivative is too complicated. I agree that the easiest approach is to specify the Stokes drift as $u^S(x,y,z,t)$ and compute the required values and derivatives internally. The difference between the Stokes drift at a specific depth and a depth-averaged Stokes drift around that depth should be small if we have sufficient resolution.

Copy link
Member Author

@glwagner glwagner Dec 27, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe asking the users to specify every derivative is too complicated.

It's what we do currently 😂

We don't have to choose to do anything specific. We can keep the existing structure, and add a new one wherein users specify the Stokes drift only. Then the test of time shall arbitrate.

Mainly I think it could save us time and effort in the end to have just one implementation. Then we know that we are all using the same code and more likely to find bugs / benefit from improvements together.

The difference between the Stokes drift at a specific depth and a depth-averaged Stokes drift around that depth should be small if we have sufficient resolution.

I do think that either way, sufficient resolution renders either discretization similar. But it would be interesting to understand the implications of inadequate resolution too.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think as @BrodiePearson pointed out we will need to specify three more spatial derivatives ($\partial_x u^s$, $\partial_y v^s$, and $\partial_z w^s$)...

I agree. I was not saying we should ignore the differences. But I cannot think of a clean way to specify the Stokes drift. I think we need both the Stokes drift at a specific depth and the depth-averaged Stokes drift around that depth. Maybe asking users to specify both is another option? I was also thinking of an option to specify the wave spectrum, or a simplified one with less frequency bins as Brandon did in MOM6. But maybe these different options can be wrapped into another layer of interface in the StokesDrift module, or defined by users. Internally in Oceananigans, I guess we need $u^s$ and $v^s$ defined as depth-averaged values at the locations of $u$ and $v$, $w^s$ defined at the location of $w$, spatial derivatives should also be depth-averaged over a grid cell?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can compute the finite volume derivatives of the Stokes drift field in the same way that we compute the discrete finite volume derivatives of the prognostic velocity field. The simplest stencil is second-order accurate (it looks just like a finite difference stencil). You could be right though that the derivatives may be known with higher accuracy if users are allowed to supply them.

It wouldn't be too much work to put an interface where the needed derivatives may be optionally provided, and otherwise are computed with a finite volume method.

end
end

subtract_stokes_drift!(model, t=time(model)) = modulate_by_stokes_drift!(model, -1, t)
add_stokes_drift!(model, t=time(model)) = modulate_by_stokes_drift!(model, +1, t)

"""
calculate_pressure_correction!(model::NonhydrostaticModel, Δt)

Calculate the (nonhydrostatic) pressure correction associated `tendencies`, `velocities`, and step size `Δt`.
"""
function calculate_pressure_correction!(model::NonhydrostaticModel, Δt)

subtract_stokes_drift!(model)

# Mask immersed velocities
foreach(mask_immersed_field!, model.velocities)

fill_halo_regions!(model.velocities, model.clock, fields(model))

solve_for_pressure!(model.pressures.pNHS, model.pressure_solver, Δt, model.velocities)

fill_halo_regions!(model.pressures.pNHS)

add_stokes_drift!(model, time(model) + Δt)

return nothing
end

Expand Down
20 changes: 15 additions & 5 deletions src/StokesDrifts.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@ using Oceananigans.Operators
using Oceananigans.Grids: AbstractGrid, node
using Oceananigans.Utils: prettysummary

using KernelAbstractions: @kernel, @index

import Adapt: adapt_structure

#####
Expand Down Expand Up @@ -148,7 +150,10 @@ const c = Center()
@inline z_curl_Uˢ_cross_U(i, j, k, grid, sw::USDnoP, U, time) = (- ℑxzᶜᵃᶠ(i, j, k, grid, U.u) * sw.∂z_uˢ(znode(k, grid, f), time)
- ℑyzᵃᶜᶠ(i, j, k, grid, U.v) * sw.∂z_vˢ(znode(k, grid, f), time))

struct StokesDrift{P, VX, WX, UY, WY, UZ, VZ, UT, VT, WT}
struct StokesDrift{P, US, VS, WS, VX, WX, UY, WY, UZ, VZ, UT, VT, WT}
uˢ :: US
vˢ :: VS
wˢ :: WS
∂x_vˢ :: VX
∂x_wˢ :: WX
∂y_uˢ :: UY
Expand All @@ -161,7 +166,10 @@ struct StokesDrift{P, VX, WX, UY, WY, UZ, VZ, UT, VT, WT}
parameters :: P
end

adapt_structure(to, sd::StokesDrift) = StokesDrift(adapt(to, sd.∂x_vˢ),
adapt_structure(to, sd::StokesDrift) = StokesDrift(adapt(to, sd.uˢ),
adapt(to, sd.vˢ),
adapt(to, sd.wˢ),
adapt(to, sd.∂x_vˢ),
adapt(to, sd.∂x_wˢ),
adapt(to, sd.∂y_uˢ),
adapt(to, sd.∂y_wˢ),
Expand Down Expand Up @@ -283,7 +291,10 @@ StokesDrift{Nothing}:
└── ∂t_wˢ: ∂t_wˢ
```
"""
function StokesDrift(; ∂x_vˢ = zerofunction,
function StokesDrift(; uˢ = zerofunction,
vˢ = zerofunction,
wˢ = zerofunction,
∂x_vˢ = zerofunction,
∂x_wˢ = zerofunction,
∂y_uˢ = zerofunction,
∂y_wˢ = zerofunction,
Expand All @@ -294,7 +305,7 @@ function StokesDrift(; ∂x_vˢ = zerofunction,
∂t_wˢ = zerofunction,
parameters = nothing)

return StokesDrift(∂x_vˢ, ∂x_wˢ, ∂y_uˢ, ∂y_wˢ, ∂z_uˢ, ∂z_vˢ, ∂t_uˢ, ∂t_vˢ, ∂t_wˢ, parameters)
return StokesDrift(uˢ, vˢ, wˢ, ∂x_vˢ, ∂x_wˢ, ∂y_uˢ, ∂y_wˢ, ∂z_uˢ, ∂z_vˢ, ∂t_uˢ, ∂t_vˢ, ∂t_wˢ, parameters)
end

const SD = StokesDrift
Expand Down Expand Up @@ -325,7 +336,6 @@ const SDnoP = StokesDrift{<:Nothing}
return wᶠᶜᶜ * (∂z_uˢ - ∂x_wˢ) - vᶠᶜᶜ * (∂x_vˢ - ∂y_uˢ)
end


@inline function y_curl_Uˢ_cross_U(i, j, k, grid, sw::SD, U, time)
wᶜᶠᶜ = ℑyzᵃᶠᶜ(i, j, k, grid, U.w)
uᶜᶠᶜ = ℑxyᶜᶠᵃ(i, j, k, grid, U.u)
Expand Down
Loading