`CFL` diagnostic is not very useful, and we don't need to depend on Tullio

[glwagner]

The CFL object:

https://github.com/CliMA/Oceananigans.jl/blob/38ea3e4c63e1a0638362b597aff4bc09eb811cdc/src/Diagnostics/cfl.jl#L9-L12

isn't really useful (there seems to be little reason to preconstruct this object). A function might be more useful, eg

cfl(simulation) = simulation.Δt / cell_advection_timescale(simulation.model)

or

advective_cfl(simulation) = ...
diffusive_cfl(simulation) = ...

Relatedly,

https://github.com/CliMA/Oceananigans.jl/blob/2d6ccfe94ce2c68857f70b2f8839020049932e00/src/Utils/cell_advection_timescale.jl#L4-L15

is a hack that's not accurate (overly restrictive) on stretched grids. I think we can replace this with

using Oceananigans.AbstractOperations: Δx, Δy, Δz

function cell_advection_timescale(grid, u, v, w)
    arch = architecture(u)

    max_u_Δx = maximum(abs, u / Δx)
    max_v_Δy = maximum(abs, v / Δy)
    max_w_Δz = maximum(abs, w / Δz)

    return 1 / max(max_u_Δx, max_v_Δy, max_w_Δz)
end

Then we can delete the infamous "accurate_cell_advection_timescale":

https://github.com/CliMA/Oceananigans.jl/blob/2d6ccfe94ce2c68857f70b2f8839020049932e00/src/Diagnostics/cfl.jl#L83-L102

For reference, the new implementation of cell_advection_timescale uses GridMetricOperation via

julia> using Oceananigans.AbstractOperations: Δx
[ Info: Precompiling Oceananigans [9e8cae18-63c1-5223-a75c-80ca9d6e9a09]
[ Info: Oceananigans will use 8 threads

julia> using Oceananigans

julia> grid = RegularRectilinearGrid(size=(2, 2, 2), extent=(1, 1, 1))
RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
        size (Nx, Ny, Nz): (2, 2, 2)
        halo (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.5, 0.5, 0.5)

julia> u = XFaceField(grid)
Field located at (Face, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (2, 2, 2)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=2, Ny=2, Nz=2)
└── boundary conditions: west=Periodic, east=Periodic, south=Periodic, north=Periodic, bottom=ZeroFlux, top=ZeroFlux, immersed=ZeroFlux

julia> u / Δx
BinaryOperation at (Face, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=2, Ny=2, Nz=2)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
└── tree: 
    / at (Face, Center, Center)
    ├── Field located at (Face, Center, Center)
    └── Δxᶠᶜᵃ at (Face, Center, Center)

@navidcy and @tomchor might have something to add about this change to the user interface.

[tomchor]

I never understood why a CFL object had to be pre-constructed, but I always assumed there was a good (probably GPU-related) reason for that. If there isn't, then I completely agree it's much simpler (without, I think, loss of functionality) to have it be two simple functions (advective and diffusive cfl).

is a hack that's not accurate (overly restrictive) on stretched grids. I think we can replace this with

using Oceananigans.AbstractOperations: Δx, Δy, Δz

function cell_advection_timescale(grid, u, v, w)
    arch = architecture(u)

    max_u_Δx = maximum(abs, u / Δx)
    max_v_Δy = maximum(abs, v / Δy)
    max_w_Δz = maximum(abs, w / Δz)

    return 1 / max(max_u_Δx, max_v_Δy, max_w_Δz)
end

I thought the reason we hadn't implemented the above is because it's slow to compute, no?

[glwagner]

If it's slow, that's because of some implementation issue rather than an intrinsic reason. We should be able to do fast reductions of abstract operations.

[glwagner]

I never understood why a CFL object had to be pre-constructed, but I always assumed there was a good (probably GPU-related) reason for that.

I guess it's helpful if we have competing functions for computing time-scales (for performance reasons). But if we have just one fast function then it doesn't seem like there's a good reason.

[tomchor]

If it's slow, that's because of some implementation issue rather than an intrinsic reason. We should be able to do fast reductions of abstract operations.

Then I see no reason not to pursue the proposed changes :)

On the user interface side: I believe "diffusive CFL" isn't quite correct, since the Courant-Friedrichs-Lewy condition is the advective one. (Although I agree that "diffusive CFL" is pretty easy to understand and intuitive.) Should we try to come up with a different name?

[glwagner]

Some of the lessons over on https://github.com/CliMA/Oceananigans.jl/issues/2024 could warrant using an object after all. The reason is a bit technical... but because of how the compiler works, we have

julia> using Oceananigans.AbstractOperations: Δx

julia> using Oceananigans

julia> grid = RegularRectilinearGrid(size=(2, 2, 2), extent=(1, 1, 1));

julia> u = XFaceField(grid);

julia> udx1 = u / Δx
BinaryOperation at (Face, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=2, Ny=2, Nz=2)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
└── tree: 
    / at (Face, Center, Center)
    ├── Field located at (Face, Center, Center)
    └── Δxᶠᶜᵃ at (Face, Center, Center)

julia> udx2 = u / Δx
BinaryOperation at (Face, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=2, Ny=2, Nz=2)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
└── tree: 
    / at (Face, Center, Center)
    ├── Field located at (Face, Center, Center)
    └── Δxᶠᶜᵃ at (Face, Center, Center)

julia> udx1 === udx2
false

The reason we get different types is because of the different "identity" functions introduced in

https://github.com/CliMA/Oceananigans.jl/blob/6dcf6ffd1fdb2febaf0d20dfab85ef1d3f83d758/src/Operators/interpolation_utils.jl#L29-L66

which we do to avoid the compiler from detecting recursion during compilation.

A few things to digest here but basically we want to cache the operations u / Δx, etc because reforming them for every calculation may incur high compilation costs.