CliMA / Oceananigans.jl

🌊 Julia software for fast, friendly, flexible, ocean-flavored fluid dynamics on CPUs and GPUs
https://clima.github.io/OceananigansDocumentation/stable
MIT License
999 stars 196 forks source link

Interface to build and output boundary conditions #3774

Open josuemtzmo opened 2 months ago

josuemtzmo commented 2 months ago

I think it will be extremely helpful to be able to output boundary conditions and other fields (like fluxes, evaporation, and others) at the ocean interface. Following the discussion from #3081, something on those lines could be implemented, but it ideally it will be nice to have something like model.boundary_conditions.

using Oceananigans.BoundaryConditions: getbc
using Oceananigans: fields

# Boundary condition extractor in "kernel function form"
@inline kernel_getbc(i, j, k, grid, boundary_condition, clock, fields) =
    getbc(boundary_condition, i, j, grid, clock, fields)

# Kernel arguments
grid = model.grid
clock = model.clock
model_fields = merge(fields(model), model.auxiliary_fields)
u, v, w = model.velocities
u_bc = u.boundary_conditions.bottom
v_bc = v.boundary_conditions.bottom

# Build operations
u_bc_op = KernelFunctionOperation{Face, Center, Nothing}(kernel_getbc, grid, u_bc, clock, model_fields)
v_bc_op = KernelFunctionOperation{Center, Face, Nothing}(kernel_getbc, grid, v_bc, clock, model_fields)

However it will be great to have something more general to simplify diagnosing simulations, that outputs the relevant boundary condition, in the example will be the bottom boundary.

josuemtzmo commented 2 months ago

I'm not sure where to start to implement this, so if you have ideas, it will be really helpful.

glwagner commented 2 months ago

Let's start with syntax which will determine where the code will live. I think it does have to be Models in order to support function boundary conditions. So something like

τx = boundary_condition(model, :u, :top)

(assuming / implying that the top bc is a flux).

We could also call this boundary_condition_operation. However this would be a misnomer in the case that the boundary condition is in fact a Field itself. Maybe not a big deal?

As for implementation, we should use dispatch to handle the various cases: we only need a KernelFunctionOperation for function boundary conditions. Otherwise the field constructor may suffice to handle constants, arrays, and fields. We would also need something for FieldTimeSeries.

But let's agree on syntax first then we can move on to those details!

josuemtzmo commented 2 months ago

Another name could be boundary_condition_field, since despite the fact that generally we will like a function to generate the boundary conditions, the return value will be a field to output. Right?

If I understand properly, the KernelFunctionOperation will handle building the output of 1D, 2D, or 3D fields right each time the function is called, while the FieldTimeSeries will build statistics and output for a time-series right?