New abstractions for buoyancy and rotation

[glwagner]

This PR rewrites the user interface for specifying rotation rates and buoyancy parameters, such as gravitational acceleration, and the relationship between tracers and density via an equation of state.

The main purpose of this PR is to move towards a more logical, functionality-driven encapsulation of parameters.

For rotation, this PR defines a new type called AbstractRotation whose subtypes are intended to describe parameters related to the background rotation rate of the model. At the moment there is only one such type: VerticalRotationAxis, or FPlane. In the future we expect to have BetaPlane and TiltedRotationAxis or FullCoriolis. The AbstractRotation types are uniquely associated with a function that computes the x, y, and z components of "f X U", associated with the background vorticity f.

For density, equations of state, and buoyancy, a new abstract type AbstractBuoyancy is defined. This type is intended to encapsulate both a gravitational acceleration parameter (when applicable), and an equation of state that contains parameters associated with the relationship between tracer values and density.

The advantage of this design is that parameters logically associated with certain terms in the equation of motion are now associated with a particular abstract type. This permits an easy extension of the rotational of buoyancy/equation of state models currently contained in the code; for example in the future we hope to support a TracerBuoyancy type for the case that buoyancy is itself a tracer and gravitational acceleration is no longer a parameter, and perhaps a TemperatureBuoyancy type for the case that only temperature contributes to buoyancy, rather than both temperature and salinity as is currently the case.

In addition, our hydrostatic and non-hydrostatic pressure fields now have the same units and can be directly summed to obtain the total pressure.

Resolves #191.
Resolves #217.
Resolves #403.

[codecov]

Codecov Report

Merging #412 into master will increase coverage by 1.57%.
The diff coverage is 69.56%.

@@            Coverage Diff             @@
##           master     #412      +/-   ##
==========================================
+ Coverage   73.28%   74.85%   +1.57%     
==========================================
  Files          25       25              
  Lines        1426     1416      -10     
==========================================
+ Hits         1045     1060      +15     
+ Misses        381      356      -25

Impacted Files	Coverage Δ
src/Operators/Operators.jl	100% <ø> (ø)
src/output_writers.jl	80.11% <ø> (ø)	⬆️
src/TurbulenceClosures/turbulence_closure_utils.jl	62.5% <ø> (ø)
...rc/TurbulenceClosures/velocity_tracer_gradients.jl	69.72% <ø> (ø)
src/TurbulenceClosures/TurbulenceClosures.jl	100% <ø> (ø)
...rbulenceClosures/turbulence_closure_diagnostics.jl	25% <ø> (ø)
src/TurbulenceClosures/closure_operators.jl	64.44% <0%> (ø)
src/clock.jl	100% <100%> (ø)	⬆️
src/Oceananigans.jl	75% <100%> (ø)	⬆️
src/time_steppers.jl	74.12% <100%> (-0.18%)	⬇️
... and 11 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a9e34a7...7fe8a47. Read the comment docs.

[ali-ramadhan]

This PR looks great to me! Should especially help when we want to move away from the f-plane (#28).

Might be good to discuss in person with @jm-c and @christophernhill to get their thoughts on this PR. I'll keep from explicitly approving this PR for now.

Kind of surprised tests passed this easily! Also nice to see new features resulting in code size going down!

[ali-ramadhan]

Hmmm, do we really need extra operators named specifically for pressure? Seems fine to me to call δx_c2f(grid, p, i, j, k) / grid.Δx, etc.

But maybe this point is moot as the operators will all be replaced by PR #283 at some point.

[ali-ramadhan]

Maybe just ∂x etc. would be a better name as it's not particular to just pressure.

[glwagner]

The operator for pressure serves no functional purpose, as you point out --- it is purely introduced to enhance readability. We certainly can use ∂x_faa (which is defined in turbulence_closures/closure_operators.jl). I think we should use functions as much as possible, so I'd prefer using a function here over hard-coding the derivative --- if this pattern is consistent with the finite volume implementation of these equations.

In TurbulenceClosures, I introduced a lot of such functions (see ∂x_u, etc) in the hope that it would enhance the readability and clarity of the hyper-complicated functions involved in the turbulence closures.

I adopted the same philosophy here.

Some improvement could be had by defining ∂x_p as an alias via

const ∂x_p = ∂x_faa

To avoid a proliferation of redundant functions (I think I did this in TurbulenceClosures, but forget when I was implementing it here).

[ali-ramadhan]

Maybe it's a GFD thing but would FPlane be a better default name than VerticalRotationAxis? I find it's pretty common to say f-plane but not very common to say vertical rotation axis to describe a constant f scenario.

[glwagner]

That's fine with me. I started out that way, but decided that VerticalRotationAxis was more descriptive and precise and thus deserved to fill the primary role here (similar to using Value instead of Dirichlet). But if you would like to change back to FPlane, go ahead.

[ali-ramadhan]

Hmmm, it's somewhat concerning that we're starting to mix turbulence closure operators with regular operators (when the two should really be merged), but this is probably fine for now and will be resolved by PR #283 when we get to it.

[glwagner]

There is no interpolation operator for function within the "regular" operators, however... ?

[ali-ramadhan]

🚀

[glwagner]

We've decided:

• No 'VerticalRotationAxis' --- should be FPlane
• model.coriolis rather than model.rotation?
• g -> gravitational_acceleration
• Implement the nonlinear equation of state proposed by Roquet et al 2015 to show an example of a non-'linear' equation of state
• require AbstractBuoyancy types to define a function called Nsq or vertical_buoyancy_gradient

[glwagner]

Notes:

• buoyancy_perturbation is sufficient for nonlinear equations of state
• Equations of state must define thermal_expansion and haline_contraction functions
• buoyancy gradient / N^2 should be calculated using these functions, in general
• we do not need total_buoyancy
• nonlinear equations of state will also be associated with a 'density anomaly' function, and a reference density.

[ali-ramadhan]

Cool stuff!

Do you think this will be ready to merge soon-ish? From implementing the nonlinear equation of state, seems that the features/API needed has become clearer.

We should probably add tests for the Roquet and others (2015) equations of state. Not sure what to test exactly though. At the very least, they should run and not crash. Beyond that could compare with hand calculations but that seems kind of useless?

PS: Sorry I've committed a great sin and edited your original comment to include "Resolves #191."

[ali-ramadhan]

Does it matter that we've been using ρ₀=1027? I guess not as this only applies to RoquetIdealizedNonlinearEOS while we've been using 1027 with LinearEquationOfState.

[glwagner]

ρ₀ has no affect on dynamics unless you're using a nonlinear equation of state (or if you use it for thermodynamics calculations).

We will need to put some thought in how to ensure that, if users do use ρ₀ for boundary conditions or thermodynamics calculations, that the same value is used for both the equation of state and those thermodynamics calculations (which in principle could occur even before declaring Model).

Possibly we can write up some 'Density' type that takes an equation of state as an argument and knows how to do the right thing (requires ρ₀ when given a LinearEquationOfState, or uses the default if its not specified and given a nonlinear equation of state.)

[ali-ramadhan]

Hmmm, yeah seems like the proper thing to do is to construct the buoyancy struct before you need to use the value of ρ₀ and then pass buoyancy to the Model constructor.

It's not the prettiest solution though.

[glwagner]

Do you think this will be ready to merge soon-ish?

Yes I think so.

The main clarifications here are that

1. The concept of 'buoyancy' is on solid ground for any equation of state, provided we are within the Boussinesq approximation

2. The dynamic component of the buoyancy gradient (the squared buoyancy frequency) should be calculated using the local thermal expansion coefficient and haline contraction coefficient (the analytical derivatives of the equation of state with respect to temperature and salinity at constant depth) --- This is (apparently) a more accurate approximation to the buoyancy gradient than calculating the buoyancy gradient with finite differences of buoyancy at different vertical levels. But this means that all equations of state should define functions that return these coefficients. In addition, we will supply a function to calculated the buoyancy frequency squared (I'm currently calling it N2).

3. For Boussinesq equations of state, we do not require the non-static components of pressure as input; instead we use the geopotential depth (corresponding to the use of hydrostatic pressure, ρ₀ g depth, to calculate compressive / thermobaric effects on density). At the moment we thus only require -z as an input --- though we may need to generalize the buoyancy implementation in the future if for some reason calculating the 'geopoential depth' requires more information (like free surface elevation). Luckily, this is a small change that just requires changing the function arguments to buoyancy_perturbation; etc.

4. The name buoyancy_perturbation is ok. It is actually potentially confusing, since "buoyancy" in ordinary usage is already associated explicitly with density anomalies --- rather than static, typically constant or z-dependent components of the density profile. So I think calling this function simply buoyancy might be better. On the other hand, we probably want to use model.buoyancy to refer to the model field that holds buoyancy-related parameters. Because of that we may want to stick with buoyancy_perturbation. Any ideas / alternative suggestions? An alternative is to perform a translation in the kernel functions; aka use buoyancy_params in the function signature of kernel(s) that call the function buoyancy.

[glwagner]

What are feelings about the nomenclature "IdealizedNonlinearEOS" versus "IdealizedNonlinearEquationOfState"?

[glwagner]

Gregg, Ali, can we get variable delz in at the same time? John

I'm happy to start working on that in October. I think we should reverse the k-index first, because I think at least part of the variable delz algorithm will depend on the direction of k.

@ali-ramadhan on buoyancy versus buoyancy_perturbation --- since this is an internal issue (not part of the API), I'm fine to stick with buoyancy_perturbation for now and continue to discuss what the best nomenclature is after this PR is merged.

[ali-ramadhan]

The point I was trying to discuss was the fact that the usual definition of "buoyancy" typically isolates the dynamic component of the force associated with density. ρg is the force density associated with gravity (buoyancy could be a force, but is typically defined with units of acceleration). But ρg does not induce motion when ρ is constant (for example), because in that case it is balanced by hydrostatic pressure (it does not contribute to dynamics).

Are you saying that you disagree and you believe the usual definition of buoyancy is ρg? Perhaps this is true --- I just want understand what your point is here.

This is probably just me being ignorant or maybe just different conventions in different fielss but I guess I remember seeing ρg being introduced as the buoyant force (per unit volume) in fluid mechanics 101. So to me, I thought this is the definition of buoyancy. But it sounds like when people say "the buoyancy" they are referring to the buoyant acceleration.

This is the key line...

To me it seems that buoyancy_perturbation(i, j, k, grid, buoyancy, C) and buoyancy_perturbation(i, j, k, grid, buoyancy_parameters, C) are somewhat similar as you can "unpack" buoyancy_parameters = model.buoyancy.

[ali-ramadhan]

on buoyancy versus buoyancy_perturbation --- since this is an internal issue (not part of the API), I'm fine to stick with buoyancy_perturbation for now and continue to discuss what the best nomenclature is after this PR is merged.

Completely agree. @suyashbire1 and I are hoping this gets merged ASAP because all our scripts are depending on old versions right now (which have other issues) and we need to update them but can't until the buoyancy and rotation API is stable (well we can, but we want to refactor all our scripts at once).

[ali-ramadhan]

a few issues to resolve before merging

Might also add changing VerticalRotationAxis to FPlane.

[glwagner]

Hmm, well ρg is literally the force density exerted by the weight of the fluid. On Wikipedia, 'buoyancy' is the force induced when you displace fluid; therefore the buoyancy of a solid object is equal to

gravitational acceleration x (volume of object x density of displaced fluid - mass of the object)

which gives the forcing acting on the body (positive upwards due to the sign of the terms).

Here we deal in densities, so the buoyant forced defined in the same manner would essentially be

gravitational acceleration x ( reference density - density of the buoyant fluid )

which of course is precisely -ρ' g, since ρ' = ρ - ρ0.

I think that from its usage the meaning of buoyancy_perturbation is clear.

Something to keep in mind as we continue to discuss this topic outside this PR is that atmospheric scientists go even a step further and define 'buoyancy' as -(ρ' - <ρ'>) g / -ρ0, where <> is a horizontal average --- thus omitting the z-dependent, horizontally uniform part of the 'total' buoyancy that does not contribute to dynamics (because it is balanced by a static component of pressure). Thus using buoyancy_perturbation to denote a quantity which, for them, comprises even non-dynamic, irrelevant components of -ρ' g / -ρ0 could seem a bit odd.

[glwagner]

@ali-ramadhan please merge if tests pass.

[ali-ramadhan]

Tests are still failing but I'll fix them.

[ali-ramadhan]

Just some suggested changes but nothing important. We can clean things up later. I think it's more important to merge this and tag a new stable release.

[ali-ramadhan]

Needs another section header here maybe.

[ali-ramadhan]

I thought we weren't using types in function signatures unless it's for dispatch? (Following the YASGuide)

[glwagner]

Hmm... I think this was originally needed for dispatch because f came with a default; however it looks like it may not be needed anymore.