Adds an f-plane implementation with arbitrary direction of rotation #1892
Conversation
|
From #1372:
Yeah I agree this might make it a bit confusing. For now I didn't delete anything and created |
|
On a refactoring note, currently Which is after I also think would be useful to use something like So I was thinking of moving this to Since this would be refactoring code, I'll wait for some feedback before doing these modifications. |
|
Thanks for creating this @tomchor , I think this is a neat idea. To help me think about how this should look, could you help me find an example you want want to do this? I can imagine that maybe you would want a general fplane and general nontraditional fplane as well. I guess it depends on the physical set up. |
Sure thing! The example I'm developing this for specifically is a tilted bottom boundary layer. (#1498) Over there I'm tilting the domain by tilting the gravity vector, really. However, to "properly" tilt the domain I'd need to tilt Coriolis too. In this case I'd need a x-component to const θ_rad = 0.2 # radians
const g̃ = (sin(θ_rad), 0, cos(θ_rad)) # domain tilt
buoyancy = Buoyancy(model=BuoyancyTracer(), vertical_unit_vector=g̃) # Tilt gravity vector
coriolis = GeneralFPlane(coriolis_frequency=1e-4, rotation_axes=g̃) # Tilt rotation vector |
I'm actually not sure. All that |
This should definitely replace |
|
As of right now, as far as I can tell, The interface I implemented is a bit simpler than
Questions:
|
There was a problem hiding this comment.
Lots of nice work here @tomchor, thanks for putting this together!
|
Thoughts on your questions:
|
I've been asked to review --- is this still the state of this PR? |
No it's not. Thanks, I forgot to change that. It's now ready for review with the caveat that I'd like to wait until I get feedback on how it's working to change the docs if that's okay, as per this comment. I just wanna make sure you're all okay with the changes in functionality. But lmk if that's not okay and I'll change the docs too. |
For API this is what I suggest: Three "modes":
The code might look something like zero_if_nothing(f) = f
zero_if_nothing(::Nothing) = 0
function ConstantBackgroundRotation(FT=Float64; fx=nothing, fy=nothing, fz=nothing, rotation_rate=Ω_Earth, rotation_axis=nothing, latitude=nothing)
if latitude !=nothing
isnothing(rotation_axis) && throw(ArgumentError("Cannot specify latitude and rotation axis."))
all(isnothing.((fx, fy, fz)) || throw(ArgumentError("Cannot specify latitude and (fx, fy, fz)."))
# calculate rotation axis
end
if rotation_axis != nothing
all(isnothing.((fx, fy, fz)) || throw(ArgumentError("Cannot specify rotation_axis and (fx, fy, fz)."))
# calculate fx, fy, fz
end
fx, fy, fz = zero_if_nothing.((fx, fy, fz)) # set default fx, fy, fz
return ConstantBackgroundRotation(FT(fx), FT(fy), FT(fz))
endThere's also the possibility of a somewhat minor optimization by keeping the possibility that I think |
I'm talking about the hydrostatic approximation in a spherical coordinate system, where do not use the Cartesian components. |
Ah, I see. I'm not sure how that would factor in since I'm not familiar with the equations. But I'll trust that for now it's back to a separate constructor as we've been talking about. |
|
If all the tests pass then this PR is ready for its final review. Notes:
|
|
Ah sorry I really do think it would be confusing to users that one has to use Sorry this is becoming laborious... if the name is changed to include Cartesian then we can merge this and discuss further in an issue. I think we actually need two types (one for nonhydrostatic and one for hydrostatic) along with a design that generalizes to any coordinate system (storing rotation axis and either |
|
This is a bit complicated but it is better to have discussion now to plan what we think is best moving forward. It does seem that we want to make sure that we can't use f- and beta-planes in a spherical geometry, and don't want to use the spherical Coriolis force in a Cartesian geometry. This can be done in a variety of ways I suppose. If adding Cartesian in the name helps for the moment then why not? That way we can merge sooner rather than later and then discuss this as an issue. We all agree that we need different Coriolis functions for Cartesian and Spherical coordinates. Why do we need different functions or nonhydrostatic (and shallow water I presumme) and hydrostatic? |
Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>
I do agree it's a bit confusing, but IMO it's also obscure to name it
I'm okay with this. It seems like the background rotation implementation could use a big refactoring, which is kind of outside the scope of this PR anyway. |
I'm also not clear, although I'm assuming that it has to do with |
Yes, when we make the hydrostatic approximation, we assume that the aspect ratio is thin and H/L is small. The hydrostatic assumption specifically refers to the use of this scaling in the vertical momentum equation (reducing it to a diagnostic equation for hydrostatic pressure). This same scaling applied to the Coriolis force leads to the "traditional" approximation such that Coriolis terms involving the vertical velocity are neglected from the horizontal momentum equations (likewise, the terms involving the horizontal momentum are neglected from the vertical momentum balance; neglecting these terms must be made consistently for the system to conserve kinetic energy). This thin-aspect-ratio approximation (probably best to avoid implicating "tradition" in model formulation...) also needs to be invoked to justify "f-plane" and "beta-plane" approximations to the Coriolis term when the numerical model is supposed to approximate oceanic motion away from the poles. (The so-called "non-traditional" terms --- the projection of the planetary vorticity into horizontal directions --- have been variously found to have a small effect on turbulent boundary layer motions. This is probably because the effect of Coriolis forces is most pronounced at the largest scales, and the effect of the horizontal Coriolis components on large scale vertical motions is suppressed by the presence of an impenetrable surface at the top and bottom.) In Oceananigans, we have to express this notion with a type to avoid adding spurious terms to the horizontal momentum equations that depend on the vertical velocity. |
|
@glwagner : I understand what you are saying about the traditional approximation. If we include the non-tradidtional terms, we obtain the Quasi-Hydrostatic model, as discussed in the MITgcm. This is something we can solve in the I'm not suggesting we change the name of the model but I thought we could use this in either model, but maybe there's a problem because of the pressure solve, which might not allow for this? |
We could indeed implement such a model; it would not be very difficult. It would mean that the hydrostatic pressure depends on the Such a modification to the hydrostatic pressure integral is also needed to introduce surface wave effects via the Craik-Leibovich approximation to a hydrostatic model (this effects can be interpreted as a modification to the background rotation rate of the fluid, with the vertical derivative of the Stokes drift acting as the Coriolis vector). |
Okay, so a longer term plan. For the moment we need that |
Note too that the problem is really specific to a spherical geometry. In a rectangular geometry, a thin aspect ratio approximation is mathematically equivalent to a change in the axis of rotation. In a spherical geometry this is no longer possible and we need to explicitly state that we are making a hydrostatic approximation. |
|
Added some dynamics tests for Coriolis with this last commit. It two a 0-D case for half an inertial period with a rotation about the There's one part that tests if the total velocity magnitude is approximately unchanged (magnitude=1), which relies on an implicit arbitrary tolerance which might be bad. I'd curious about your feedback on that one.
Per the comment above I'm going to change the name to |
For physics tests we can't avoid introducing an arbitrary tolerance. So its ok. That's one reason why physics tests in CI are a bit problematic and we also need validation tests analyzed by humans. |
Np. This was a good illustration of scientific debate haha |
Closes #1372
For now I just wanted to save up the calculations. Once we decide how to implement the feature I'll finish the PR.