FieldTimeSeries for time-dependent boundary condition

[simone-silvestri]

The intention of this PR is to extend the capability of FieldTimesSeries such that it can be used in time-dependent boundary conditions and forcing.

This PR develops over 2 main points, but it is still a work in progress so suggestions are welcome

1. Allow writing down files with set!ing the OnDisk flavor of a field time series so that it aligns with the format of our jld2 output writer. In this way, it is possible to easily format BC from different file types (such as binary or text) to be used with Oceananigans.
   An example of this is:

grid = LatitudeLongitudeGrid(size = (16, 16, 10), latitude = (-60, -40), longitude = (-10, 10), z = (0, 1))
u = XFaceField(grid)
fts = FieldTimeSeries{location(u)...}(grid, 1:100, backend = OnDisk(), path = "testfile.jld2", name = "u")
for i in 1:length(fts.times)
    set!(u, 2i)
    set!(fts, u, i)
end

This will generate a file called testfile.jld2 with the following structure

julia> f = jldopen("testfile.jld2")
JLDFile /Users/simonesilvestri/development/Oceananigans.jl/testfile.jld2 (read-only)
 ├─📂 serialized
 │  └─🔢 grid
 └─📂 timeseries
    ├─📂 u
    │  ├─🔢 1
    │  ├─📂 serialized
    │  │  ├─🔢 location
    │  │  └─ ⋯ (2 more entries)
    │  └─ ⋯ (99 more entries)
    └─📂 t (100 entries)

which can be easily read by the other field time series types.

2. To do linear interpolation we need at least 2 fields in memory, so OnDisk will not do. On the other hand, we might not want all fields in memory as if we are dealing with forcings that might overwhelm the memory (especially on the GPU). So the proposal is to implement a Chunked abstraction that only keeps in memory a "chunk" of the data. The details of this implementation are still open do be decided, especially if we want an automatic update of the chunk if we index into an index not existing in memory or if we want the user to be responsible in updating the data in memory through something like a callback

julia> fts3 = FieldTimeSeries("testfile.jld2", "u", backend = Chunked(chunk_size = 4))
17×16×10×100 FieldTimeSeries{Chunked} located at (Face, Center, Center) on CPU
├── grid: 16×16×10 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── indices: (Colon(), Colon(), Colon())
└── data: 23×22×16×4 OffsetArray(::Array{Float64, 4}, -2:20, -2:19, -2:13, 1:4) with eltype Float64 with indices -2:20×-2:19×-2:13×1:4
    └── max=8.0, min=0.0, mean=1.67984

This is not final and all open to suggestions/changes/improvement

maybe interesting for @yuchenma23

[glwagner]

How would we make this work with time::DateTime? I don't think we need that feature right now (it matters only if you have observational time-series that span more than one year)

[simone-silvestri]

Hmm yeah

[simone-silvestri]

Hmmm yeah, do AbstractTimes support operations?

[glwagner]

Suggested change

	t₁, t₂ = index_binary_search(fts.times, time, Ntimes)
	n₁, n₂ = index_binary_search(fts.times, time, Ntimes)

[glwagner]

most places we use n for the time index

[glwagner]

hmm maybe we can extend index_binary_search to work for x::AbstactRange. Might make sense to call it "search_index_bracket" after that or something.

[simone-silvestri]

at the moment for chunked arrays the user can change with set!. For the above example of fts3 which has a chunk size of 4:

julia> fts3[1]
17×16×10 Field{Face, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 16×16×10 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── boundary conditions: Nothing
└── data: 23×22×16 OffsetArray(view(::Array{Float64, 4}, :, :, :, 1), -2:20, -2:19, -2:13) with eltype Float64 with indices -2:20×-2:19×-2:13
    └── max=2.0, min=2.0, mean=2.0

julia> fts3[7]
ERROR: BoundsError: attempt to access 23×22×16×4 Array{Float64, 4} at index [1:23, 1:22, 1:16, 7]

julia> set!(fts3, 7:10)

julia> fts3[7]
17×16×10 Field{Face, Center, Center} on LatitudeLongitudeGrid on CPU
├── grid: 16×16×10 LatitudeLongitudeGrid{Float64, Bounded, Bounded, Bounded} on CPU with 3×3×3 halo and with precomputed metrics
├── boundary conditions: Nothing
└── data: 23×22×16 OffsetArray(view(::Array{Float64, 4}, :, :, :, 1), -2:20, -2:19, -2:13) with eltype Float64 with indices -2:20×-2:19×-2:13
    └── max=14.0, min=14.0, mean=14.0

julia> fts3[1]
ERROR: BoundsError: attempt to access 23×22×16×4 Array{Float64, 4} at index [1:23, 1:22, 1:16, -5]
Stacktrace:

So at the moment the chunked field does not "auto-load". I am not sure whether it is better to autoload or not.

[glwagner]

What's the point of AbstractSolver?

[simone-silvestri]

It was a test I was doing, I found a better way to solve the same problem

[glwagner]

What is remaining to do on this PR? Do you need help @simone-silvestri ?

[simone-silvestri]

Should be ready. Maybe a couple of tests

[simone-silvestri]

Should be ready to go

[glwagner]

I suggest simplifying the validation test to use FieldTimeSeries just for the top temperature boundary condition. Also, making path and name properties of FieldTimeSeries so they aren't replicated for each backend.

[simone-silvestri]

Suggestions implemented! I ll merge this PR