(0.8.8) Add DateTime clock timestepping test for OceanSeaIceModel + fix import issues with Oceananigans v0.102

[navidcy]

This PR adds a timestepping test for the OceanSeaIceModel. It also fixes a bunch of import-related errors that have to do with ExplicitImport changes in Oceananigans (I hope I did that correctly; cc @giordano).

This PR depends on JuliaRegistries/General#143185 because it needs Oceananigans v0.102.1 that includes CliMA/Oceananigans.jl#4885

[navidcy]

This PR supersedes #687

[giordano]

It also fixes a bunch of import-related errors that have to do with ExplicitImport changes in Oceananigans (I hope I did that correctly; cc @giordano).

Yes, those changes look correct. I think the Time change is a good one, it should really be in Units, not Utils, the prettysummary one is more unfortunate though, I believe it should be in Utils, not Grid, but I didn't find a way to move the modules around. However we can have the prettysummary symbol be available in Utils too, not to break downstream packages, if you think it's a good idea.

[simone-silvestri]

Please, do have a prettysummary symbol in Utils! This change is breaking all the downstream packages

[navidcy]

Seems like some GPU tests are failing surpassing the kernel 4kb limit?

[simone-silvestri]

apparently that started happening as we switched to julia 1.12

[simone-silvestri]

Let's merge this since it is blocking further development. We can figure out the GPU issue of the parameter size in a new PR. In general, I think that we should look at stopping a bit the development and fix all the tests of ClimaOcean.
I would not do it though until Oceananigans is a bit more stable.

[simone-silvestri]

Note we also need to change buoyancy to buoyancy_force

[navidcy]

The tests with PrescribedAtmosphere and DateTime clock break on GPU; something related to scalar indexing in the TimeInterpolator.

https://buildkite.com/clima/climaocean-ci/builds/5280#019afb3d-9da1-43de-aa67-0fad6d50f04d/187-604

Any idea why? In particular:

using ClimaOcean, Oceananigans, CUDA, Dates

start_time = Dates.DateTime(1993, 2, 1)
Δt_seconds = 3050.0
Δt_period = Dates.Second(3050)
steps = 5
stop_time = start_time + steps * Δt_period
times = [start_time + m * Δt_period for m in 0:steps]

arch = GPU()

grid = LatitudeLongitudeGrid(arch; size = (1, 1, 20), longitude = (10, 11), latitude = (20.2, 21.2), z = (-200, 0), topology = (Bounded, Bounded, Bounded))

ocean_clock = Clock(time=start_time)
ocean = ocean_simulation(grid; clock=ocean_clock, Δt=Δt_seconds, tracers=(:T, :S), verbose=false, closure=nothing)
atmosphere_clock = Clock(time=start_time)
atmosphere = PrescribedAtmosphere(grid, times; clock=atmosphere_clock)
radiation = Radiation(arch)
coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)

ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] errorscalar(op::String)
    @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex
    @ /g/data/v45/nc3020/.julia/packages/GPUArrays/0F4Dn/src/host/indexing.jl:50 [inlined]
  [6] searchsortedfirst(v::CuArray{DateTime, 1, CUDA.DeviceMemory}, x::DateTime, lo::Int64, hi::Int64, o::Base.Order.ForwardOrdering)
    @ Base.Sort ./sort.jl:192
  [7] searchsortedfirst
    @ ./sort.jl:308 [inlined]
  [8] searchsortedfirst
    @ ./sort.jl:310 [inlined]
  [9] find_time_index
    @ /g/data/v45/nc3020/.julia/dev/Oceananigans/src/OutputReaders/field_time_series_indexing.jl:115 [inlined]
 [10] interpolating_time_indices
    @ /g/data/v45/nc3020/.julia/dev/Oceananigans/src/OutputReaders/field_time_series_indexing.jl:77 [inlined]
 [11] TimeInterpolator
    @ /g/data/v45/nc3020/.julia/dev/Oceananigans/src/OutputReaders/field_time_series_indexing.jl:35 [inlined]
 [12] interpolate_atmosphere_state!(interfaces::ComponentInterfaces{…}, atmosphere::PrescribedAtmosphere{…}, coupled_model::OceanSeaIceModel{…})
    @ ClimaOcean.OceanSeaIceModels.InterfaceComputations /g/data/v45/nc3020/ClimaOcean.jl/src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl:81
 [13] #update_state!#4
    @ /g/data/v45/nc3020/ClimaOcean.jl/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl:46 [inlined]
 [14] update_state! (repeats 2 times)
    @ /g/data/v45/nc3020/ClimaOcean.jl/src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl:37 [inlined]
 [15] initialization_update_state!(model::OceanSeaIceModel{FreezingLimitedOceanTemperature{…}, PrescribedAtmosphere{…}, Simulation{…}, ComponentInterfaces{…}, Clock{…}, GPU{…}})
    @ ClimaOcean.OceanSeaIceModels /g/data/v45/nc3020/ClimaOcean.jl/src/OceanSeaIceModels/ocean_sea_ice_model.jl:49
 [16] OceanSeaIceModel(ocean::Simulation{…}, sea_ice::FreezingLimitedOceanTemperature{…}; atmosphere::PrescribedAtmosphere{…}, radiation::Radiation{…}, clock::Clock{…}, ocean_reference_density::Float64, ocean_heat_capacity::Float64, sea_ice_reference_density::Int64, sea_ice_heat_capacity::Int64, interfaces::Nothing)
    @ ClimaOcean.OceanSeaIceModels /g/data/v45/nc3020/ClimaOcean.jl/src/OceanSeaIceModels/ocean_sea_ice_model.jl:206
 [17] top-level scope
    @ REPL[24]:1
Some type information was truncated. Use `show(err)` to see complete types.

Any ideas why? Is the test faulty (e.g. creating times as an Array of dates) or is there an issue with the searchsortedfirst?

[glwagner]

I think that there's a part of the Oceananigans code that tries to do a computation on the CPU. Maybe dispatch is slipping up and missing it. Not sure why it wouldn't be tested in Ocewananigans. It's something like cpu_time_interpolation_indices

[glwagner]

https://github.com/CliMA/Oceananigans.jl/blob/58c6b7a7b4760abc4da76ddef1ae8d198ea82b6b/src/OutputReaders/field_time_series_indexing.jl#L115

doesn't this work on GPU?

[glwagner]

I think you can create a simpler MWE by just launching a kernel that merely tries to searchsortedfirst on a CuArray of dates. If that doesn't work it could become an issue on GPUArrays

[navidcy]

julia> using CUDA

julia> searchsortedfirst(CuArray([1, 2, 4, 5, 5, 7]), 5)
ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:44
 [2] errorscalar(op::String)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
 [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
 [4] assertscalar(op::String)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
 [5] getindex
   @ /g/data/v45/nc3020/.julia/packages/GPUArrays/0F4Dn/src/host/indexing.jl:50 [inlined]
 [6] searchsortedfirst(v::CuArray{Int64, 1, CUDA.DeviceMemory}, x::Int64, lo::Int64, hi::Int64, o::Base.Order.ForwardOrdering)
   @ Base.Sort ./sort.jl:192
 [7] searchsortedfirst
   @ ./sort.jl:308 [inlined]
 [8] searchsortedfirst(v::CuArray{Int64, 1, CUDA.DeviceMemory}, x::Int64)
   @ Base.Sort ./sort.jl:310
 [9] top-level scope
   @ REPL[107]:1

[navidcy]

Somehow cpu_time_interpolation_indices is never called..

[glwagner]

julia> using CUDA

julia> searchsortedfirst(CuArray([1, 2, 4, 5, 5, 7]), 5)
ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:44
 [2] errorscalar(op::String)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
 [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
 [4] assertscalar(op::String)
   @ GPUArraysCore /g/data/v45/nc3020/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
 [5] getindex
   @ /g/data/v45/nc3020/.julia/packages/GPUArrays/0F4Dn/src/host/indexing.jl:50 [inlined]
 [6] searchsortedfirst(v::CuArray{Int64, 1, CUDA.DeviceMemory}, x::Int64, lo::Int64, hi::Int64, o::Base.Order.ForwardOrdering)
   @ Base.Sort ./sort.jl:192
 [7] searchsortedfirst
   @ ./sort.jl:308 [inlined]
 [8] searchsortedfirst(v::CuArray{Int64, 1, CUDA.DeviceMemory}, x::Int64)
   @ Base.Sort ./sort.jl:310
 [9] top-level scope
   @ REPL[107]:1

I think you have to call this inside a kernel.

[navidcy]

Could you confirm that as far as you understand the test is not balony? Because it was generated by AI and I am not that trustful to it yet...

[glwagner]

Could you confirm that as far as you understand the test is not balony? Because it was generated by AI and I am not that trustful to it yet...

Ok sorry I didn't grasp what you were asking. You should never "trust" AI --- no more than you should trust me! Even if I wrote the test you should criticize it. All code should be reviewed the same way no matter who made it (AI or human).

I just want to point out that this is really an opinion about the user interface. Do we want code like this to work:

start_time = Dates.DateTime(1993, 2, 1)
Δt_seconds = 3050.0
Δt_period = Dates.Second(3050)
steps = 5
stop_time = start_time + steps * Δt_period
times = [start_time + m * Δt_period for m in 0:steps]

arch = GPU()
grid = LatitudeLongitudeGrid(arch; size = (1, 1, 20), longitude = (10, 11), latitude = (20.2, 21.2), z = (-200, 0), topology = (Bounded, Bounded, Bounded))

atmosphere_clock = Clock(time=start_time)
atmosphere = PrescribedAtmosphere(grid, times; clock=atmosphere_clock)

I think so.

Also does this work with Float64 times? Even more so in that case!

It can be simplified. I don't think we need to write the time-step twice (or do we? we should fix that if so)

start_time = Dates.DateTime(1992, 2, 1)
Δt = Dates.Second(3050)
Nt = 5
times = range(start_time, step=Δt, length=Nt) |> collect # convert to Array for test, even though unecessary
grid = LatitudeLongitudeGrid(GPU(); size = (1, 1, 20), longitude = (10, 11), latitude = (20.2, 21.2), z = (-200, 0))
clock = Clock(time=start_time)
atmosphere = PrescribedAtmosphere(grid, times; clock)

It is interesting that the atmosphere and ocean may want different clocks. We should not require the user to know this in my opinion.

[glwagner]

Also I think the error suggests we need a test in Oceananigans for FieldTimeSeries. I can work on that.