Skip to content

Fix some FieldTimeSeries indexing bugs #4505

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 11 commits into from
May 22, 2025
4 changes: 2 additions & 2 deletions src/OutputReaders/field_time_series.jl
Original file line number Diff line number Diff line change
Expand Up @@ -247,9 +247,9 @@ mutable struct FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, P, N, KW}
end

if times isa AbstractArray
# Try to convert to a range, cuz
# Try to convert to a lighter-weight range for efficiency
time_range = range(first(times), last(times), length=length(times))
if all(time_range .≈ times) # good enough for most
if isapprox(time_range, times)
times = time_range
end

Expand Down
23 changes: 18 additions & 5 deletions src/OutputReaders/set_field_time_series.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,23 +7,36 @@ using Oceananigans.Architectures: cpu_architecture

iterations_from_file(file) = parse.(Int, keys(file["timeseries/t"]))

find_time_index(time::Number, file_times) = findfirst(t -> t ≈ time, file_times)
find_time_index(time::AbstractTime, file_times) = findfirst(t -> t == time, file_times)
function find_time_index(time::Number, file_times, Δt)
# Accommodate round-off discrepancies between the FTS times and file times
# (see https://github.com/CliMA/Oceananigans.jl/pull/4505)
ϵ = 100 * eps(eltype(file_times))
return findfirst(t -> isapprox(t, time, atol=ϵ*Δt), file_times)
end

find_time_index(time::AbstractTime, file_times, Δt) = findfirst(t -> t == time, file_times)

function set!(fts::InMemoryFTS, path::String=fts.path, name::String=fts.name; warn_missing_data=true)
file = jldopen(path; fts.reader_kw...)
file_iterations = iterations_from_file(file)
file_times = [file["timeseries/t/$i"] for i in file_iterations]
close(file)


# Compute a timescale for comparisons
Δt = mean(diff(file_times))

arch = architecture(fts)

# TODO: a potential optimization here might be to load
# all of the data into a single array, and then transfer that
# to parent(fts).

# Index times on the CPU
cpu_times = on_architecture(CPU(), fts.times)

for n in time_indices(fts)
t = fts.times[n]
file_index = find_time_index(t, file_times)
t = cpu_times[n]
file_index = find_time_index(t, file_times, Δt)

if isnothing(file_index) # the time does not exist in the file
if warn_missing_data
Expand Down