Add a DensityField and/or PotentialDensityField similar to BuoyancyField for saving density during simulations

[jbisits]

For some experiments I am running I have created a DensityField and PotentialDensityField that saves the seawater density computed by the SeawaterBuoyancy model being used. It is computed in a similar manned to BuoyancyField in buoyancy_field.jl,

@inline function densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C)
    T, S = get_temperature_and_salinity(b, C)
    return  b.equation_of_state.reference_density + ρ′(i, j, k, grid, b.equation_of_state, T, S)
end
density(model) = density(model.buoyancy, model.grid, model.tracers)
density(b, grid, tracers) =
    KernelFunctionOperation{Center, Center, Center}(densityᶜᶜᶜ, grid, b.model, tracers)
DensityField(model) = Field(density(model))

To compute a PotentialDensityField

"Extend `ρ′` to compute at user defined reference pressure"
Oceananigans.BuoyancyModels.ρ′(i, j, k, grid, eos, θ, sᴬ, pᵣ) = ρ′(θ_and_sᴬ(i, j, k, θ, sᴬ)..., pᵣ, eos)

@inline function potential_densityᶜᶜᶜ(i, j, k, grid, b::SeawaterBuoyancy, C, parameters)
    T, S = get_temperature_and_salinity(b, C)
    pᵣ = parameters.pᵣ
    return  b.equation_of_state.reference_density + ρ′(i, j, k, grid, b.equation_of_state, T, S, pᵣ)
end
potential_density(model, parameters) = potential_density(model.buoyancy, model.grid,
                                                         model.tracers, parameters)
potential_density(b, grid, tracers, parameters) =
    KernelFunctionOperation{Center, Center, Center}(potential_densityᶜᶜᶜ, grid, b.model,
                                                    tracers, parameters)
PotentialDensityField(model, parameters) = Field(potential_density(model, parameters))

If this is something you think others may be interested in I could look to add it in a PR but if not no problem!

[glwagner]

I think something like this would be useful!

Where do these expressions come from? I think the total density is missing a compressibility term that we don't include in the density perturbation ρ′. For example, in equation 13 of Rouqet et al 2015 the variable ρ′ is called r, implying that another quantity r0 should be added to obtain the total density.

I also don't understand the motivation behind the definition of potential density (it could make sense, there are just no comments or references). I'm suspicious of the appearance of a variable called pr, which evokes "pressure": it was pointed out by Young 2010 that geopotential height rather than pressure properly enters Boussinesq expressions for density. Maybe I'm missing something.

Something more minor: in other discussions, we decided it would better encourage best practices if we define utilities that return the appropriate KernelFunctionOperation, and avoid defining utilities that define Field. This is because for most sizeable experiments the Field is undesirable (allocating a lot of memory because it's a 3D field).

I hope that isn't discouraging --- I think these objects would be useful! They should be defined in Oceananigans.Models.

[glwagner]

I think that is a great solution for computing density! We may not need the intermediate wrapper, as explained below

For potential density --- can we design an abstraction for geopotential height that would allow the user to extract potential density from the density utility?

Here's my thought:

@inline ρ(i, j, k, grid, eos, T, S, Z) = @inbounds ρ(T[i, j, k], S[i, j, k], Z[i, j, k], eos)

seawater_density(grid, eos, temperature, salinity, geopotential_height) =
    KernelFunctionOperation{Center, Center, Center}(ρ, grid, eos, temperature, salinity, geopotential_height)

const ModelsWithBuoyancy = Union{NonhydrostaticModel, HydrostaticFreeSurfaceModel}

# some nice fallbacks
model_temperature(bf, model)     = model.tracers.T
model_salinity(bf, model)        = model.tracers.S
model_geopotential_height(model) = KernelFunctionOperation{Center, Center, Center}(Zᶜᶜᶜ, model.grid)

using Oceananigans.BuoyancyModels: SeawaterBuoyancy, Zᶜᶜᶜ
using Oceananigans.Fields: field

const ConstantTemperatureSB = SeawaterBuoyancy{FT, EOS, <:Number, <:Nothing}
const ConstantSalinitySB = SeawaterBuoyancy{FT, EOS, <:Nothing, <:Number}

model_temperature(b::ConstantTemperatureSB, model) = b.model.constant_temperature
model_salinity(b::ConstantSalinitySB, model) = b.model.constant_salinity

"""
    seawater_density(model; temperature, salinity, geopotential_height)

Return a `KernelFunctionOperation` that computes the in-situ density of seawater
with (gridded) `temperature`, `salinity`, and at `geopotential_height`.
"""
function seawater_density(model::ModelsWithBuoyancy;
                          temperature = model_temperature(model.buoyancy.model, model),
                          salinity = model_salinity(model.buoyancy.model, model),
                          geopotential_height = model_geopotential_height(model))

    # Convert function or constant user input to AbstractField
    grid = model.grid
    loc = (Center, Center, Center)
    temperature = field(loc, temperature, grid)
    salinity = field(loc, salinity, grid)
    geopotential_height = field(loc, geopotential_height, grid)

    eos = model.buoyancy.model.equation_of_state
    return seawater_density(grid, eos, temperature, salinity, geopotential_height)
end

Now, users have the flexibility to construct different flavors of density by changing inputs to seawater_density. For example, to fix the geopotential height at 0 m:

σ0 = seawater_density(model, geopotential_height=0)

Any "field-like" input should work here

We should add some examples to the docstring.

What do you think @jbisits ?

[jbisits]

Looks like an excellent implementation! Much neater having a single seawater_density function to do all the hard work using the abstraction of geopotential height. In terms of test I could add one to a docstring example and some simple checks (e.g. instantiation and simple computations) in Oceananigans.jl. Given ρ lives in SewaterPolynomials.jl should any tests for the accuracy of the function ρ using reference geopotential height (against e.g. gsw_rho) be done in SeawaterPolynomials.jl? I know there are some already so this may not be necessary. If that all sounds good I can get started on this in the next day or two.

[glwagner]

should any tests for the accuracy of the function ρ using reference geopotential height (against e.g. gsw_rho) be done in SeawaterPolynomials.jl?

One could definitely calculate more "reference values" (in addition to the ones already used for tests) and I think this would be welcome. I also think a validation-type experiment could be appropriate (eg calculate global density with climatological T, S, and compare the error of the polynomial approximation to gsw_rho). All that said, I do think the minimum necessary tests are there now.

and some simple checks (e.g. instantiation and simple computations) in Oceananigans.jl

Yes, some integration test with a toy model that shows the function runs without erroring would be great! Extra credit maybe could show consistency with the SeawaterPolynomials function for a model state that's constant T, S. Or something like that...

One thing I'm not sure about is what field does when given AbstractOperation as input. That's something to think about maybe since we would want to leave abstract operations as is rather than create a field.

[jbisits]

Great, this gives me plenty to get started on. Thanks for all the help! I will open the PR when it is ready for review or if I hit any major roadblocks.

[jbisits]

Thanks @glwagner for all your help with this discussion and the implementation!