Add biogeochemical interface to ocean_simulation

[vtamsitt]

This PR addresses #107 by adding a biogeochemistry kwarg to ocean_simulation, adding user input boundary_conditions as a kwarg and moves tracer to be explicitly spelled out be the user as a kwarg.

Right now the user input boundary_conditions will be merged with hardcoded boundary conditions for T and S, unless the user supplies BCs for T and S in which case they will override the hardcoded ones. As @glwagner pointed out

it would be nice if we only replace the boundary conditions we have to, ie we keep the surface fluxes of T, S even if we need to use user-supplied T or S conditions at the bottom or on immersed boundaries

but I haven't implemented that yet.

[glwagner]

Something along these lines might work too:

tracer_advection = NamedTuple(name => tracer_advection for name in tracers)
tracer_advection = merge(tracer_advection, (; e=nothing))

[glwagner]

Also should we have something that checks if the user has already provided a named tuple? something like

if !(tracer_advection isa NamedTuple)
    # expand tracer_advection into a named tuple with `e=nothing`
end

That way people can experiment with also advecting TKE which may be useful.

cc @simone-silvestri

[glwagner]

Dictionaries also can be used which may be a little easier to understand:

tracer_advection = Dict{Symbol, Any}(name => tracer_advection for name in tracers)
tracer_advection[:e] = nothing
tracer_advection = NamedTuple(name => tracer_advection[name] for name in tracers)

[vtamsitt]

Thanks @glwagner for the suggestions

[glwagner]

🚀

[glwagner]

Looking good so far! Dictionaries can be useful for building up keyword arguments / parameters with values that have to be overwritten, since it's easy to overwrite entries in a dictionary.

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 16 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (dd9148a) to head (9cd5f7e).

Files	Patch %	Lines
src/OceanSimulations/OceanSimulations.jl	0.00%	16 Missing ⚠️

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #123   +/-   ##
=====================================
  Coverage   0.00%   0.00%           
=====================================
  Files         31      31           
  Lines       1729    1742   +13     
=====================================
- Misses      1729    1742   +13     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[vtamsitt]

@glwagner in the most recent commit I added a way to merge the user supplied boundary_conditions while keeping the hard coded T and S boundary conditions that aren't supplied by the user (ie we keep the surface fluxes of T, S even if we need to use user-supplied T or S conditions at the bottom or on immersed boundaries). It seems to work, but its ugly... any suggestions to do this more cleanly?

[glwagner]

@glwagner in the most recent commit I added a way to merge the user supplied boundary_conditions while keeping the hard coded T and S boundary conditions that aren't supplied by the user (ie we keep the surface fluxes of T, S even if we need to use user-supplied T or S conditions at the bottom or on immersed boundaries). It seems to work, but its ugly... any suggestions to do this more cleanly?

I think what you wrote captures the gist of it! But maybe we can make a few tweaks, I'll suggest something

[glwagner]

Hmm yes actually I realized we need to do something more like this:

using Oceanaingans.BoundaryConditions: DefaultBoundaryCondition

# create a utility that filters `DefaultBoundaryCondition`
boundaries = (:west, :east, :south, :north, :bottom, :top, :immersed)

"""
    field_bcs_to_dict(field_bcs)

Return a `Dict` with the properties of `field_bcs` filtered to remove `DefaultBoundaryCondition`.
"""
function field_bcs_to_dict(field_bcs)
    all_bcs = Dict(boundary => getproperty(field_bcs, boundary) for boundary in boundaries)
    non_default_bcs = filter(kv -> !isa(kv[2], DefaultBoundaryCondition), all_bcs)
    return non_default_bcs
end

# inside ocean_simulation...

# make clear these are the defaults
default_boundary_conditions = Dict(
    :u => Dict(:top => FluxBoundaryCondition(Jᵘ), :bottom => u_bot_bc),
    :v => Dict(:top => FluxBoundaryCondition(Jᵛ), :bottom => v_bot_bc),
    :T => Dict(:top => FluxBoundaryCondition(Jᵀ)),
    :S => Dict(:top => FluxBoundaryCondition(Jˢ))
)

# Build a list of the default + user-specified boundary condition names
user_bc_names = keys(boundary_conditions)
bc_names = tuple(:u, :v, :T, :S, user_bc_names...) |> unique

# Merge user and default boundary conditions
merged_boundary_conditions = Dict()

for fieldname in bc_names

    default_bcs = get(default_boundary_conditions, fieldname, Dict())

    if fieldname ∈ user_bc_names
        user_bcs = field_bcs_to_dict(boundary_conditions[fieldname])
    else
        user_bcs = Dict()
    end

    # Generate a dictionary specifying `FieldBoundaryConditions` for each `fieldname`
    bc_dict = merge(default_bcs, user_bcs)
    merged_boundary_conditions[fieldname] = FieldBoundaryConditions(; bc_dict...)
end

# Convert to NamedTuple
merged_boundary_conditions = NamedTuple(merged_boundary_conditions)

A few things we could do (in a future PR) are:

1. Print some information if the user overwrites the default boundary conditions
2. Make sure boundary_conditions is a NamedTuple and also populated correctly with FieldBoundaryConditions

[glwagner]

Not super pretty either...

[glwagner]

Realized my code sketch was missing the conversion to FieldBoundaryConditions... could be missing some other things too...

[vtamsitt]

Realized my code sketch was missing the conversion to FieldBoundaryConditions... could be missing some other things too...

Ok no worries. I'll go through and test it and make sure it's all working, thanks.

[vtamsitt]

I dropped the ball on this but picking it back up now and going to test it and hopefully it'll be ready for review this week!