Linear time interpolation in `FieldTimeSeries` #3236

simone-silvestri · 2023-08-26T15:34:46Z

Allow linear time interpolation between different indices of the field time series triggered by a Float64 (AbstractTime not supported at the moment) index.

julia> fts2 = FieldTimeSeries("testfile.jld2", "u");

julia> fts2[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> fts2[2]
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}, :, :, :, 2), -2:20, -2:19, -2:13) with eltype Float64 with indices -2:20×-2:19×-2:13
    └── max=4.0, min=4.0, mean=4.0

julia> fts2[1.35]
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: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: ZeroFlux, north: ZeroFlux, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: BinaryOperation at (Face, Center, Center)
├── status: time=0.0
└── data: 23×22×16 OffsetArray(::Array{Float64, 3}, -2:20, -2:19, -2:13) with eltype Float64 with indices -2:20×-2:19×-2:13
    └── max=2.7, min=2.7, mean=2.7

julia> fts2[1, 1, 1, 13]
26.0

julia> fts2[1, 1, 1, 14]
28.0

julia> fts2[1, 1, 1, 13.67]
27.34

navidcy · 2023-08-28T14:54:29Z

Question: If time is DateTime then only integer indices work but it fails for fractional indices. Correct?

simone-silvestri · 2023-08-28T14:58:41Z

Well, fractional indices work, just not DateTimes indices

I am not super familiar with what date-times supports in terms of operations, if it supports summation and division it is straightforward to extend indexing with to DateTimes

navidcy · 2023-08-28T15:48:58Z

I think it can be done but it’s definitely not as straightforward as +, /.

Happy to approve now and open issue to generalize the method for TimeDates later.

glwagner · 2023-08-28T16:54:55Z

Are we ok with this as a user interface? it's a little implicit. We could alternatively use a function-based API to make things a little more obvious, something like

fts[1, 2, 3, 4] # get 4th time-index

at_time(4, fts, 1, 2, 3) # linearly interpolate to t=4

mainly i'd be worried about issues like

fts[1, 2, 3, 4] \ne fts[1, 2, 3, 4.0]

which is rather easy to confuse?

We also might be able to use syntax like

fts[1, 2, 3, time=4]

if that is performant. Or

fts[1, 2, 3, Time(4)]

navidcy · 2023-08-28T17:41:37Z

Out of curiosity @simone-silvestri could you check

fts2[1, 2, 3, 4] == fts2[1, 2, 3, 4.0]

in your example with the .jld2 file above?

glwagner · 2023-08-28T18:00:39Z

Out of curiosity @simone-silvestri could you check
fts2[1, 2, 3, 4] == fts2[1, 2, 3, 4.0]
in your example with the .jld2 file above?

the outcome depends on whether fts.times = [1, 2, 3, 4] or not.

simone-silvestri · 2023-08-28T18:02:49Z

Another way would be to pass Oceananigans' Clock as an index.

This might be convenient because we usually pass the clock to functions (such as forcing and bc) so we always have it available

…to ss/rework-fts-1

simone-silvestri · 2023-08-28T18:09:02Z

@navidcy I added a way to not interpolate so that in case 4 belongs to fts2.times then fts2[4.0] == fts2[4]

but this shows a bit the confusion because 4.0 here is a time and 4 is an index.

glwagner · 2023-08-28T18:09:08Z

We wouldn't have Clock available in post-processing.

simone-silvestri · 2023-08-28T18:10:15Z

So I would propose a clock implementation for easy internal manipulation which calls the at_time function available for postprocessing

glwagner · 2023-08-28T18:10:35Z

Clock does save the definition of a new type. However, I think it is misleading because it contains extra information like the iteration and stage. Here, the time-interpolation is independent of those.

Why not a new type Time? What is your opposition to that?

glwagner · 2023-08-28T18:11:37Z

Basically I think the "convenience" of not having to define a new type is trivial. You can define Time(clock::Clock).

simone-silvestri · 2023-08-28T18:13:11Z

I am thinking about the ease of use in bc, for example

@inline time_interpolated_bc(i, j, grid, clock, fields, p) = p.time_array[i, j, 1, clock]

vs

@inline function time_interpolated_bc(i, j, grid, clock, fields, p) 
   time = Time(clock)
   return p.time_array[i, j, 1, time]
end

I find the first implementation a bit more straightforward, but we can also have both

glwagner · 2023-08-28T18:15:19Z

I am thinking about the ease of use in bc, for example
@inline time_interpolated_bc(i, j, grid, clock, fields, p) = p.time_array[i, j, 1, clock]
vs
@inline function time_interpolated_bc(i, j, grid, clock, fields, p) 
   time = Time(clock)
   return p.time_array[i, j, 1, time]
end
I find the first implementation a bit more straightforward, but we can also have both

I propose we prioritize the user interface rather than the source code implementation, which basically has to be written once then left alone

glwagner · 2023-08-28T18:17:35Z

also worth checking out

https://github.com/rafaqz/DimensionalData.jl

the Time object is a "selector" in their terminology. IF we have something similar here, we could eventually consider sub typing AbstractDimArray (which has been proposed before), which gives us a lot of powerful functionality. Someone passionate about that stuff has to take that on for a full implementation but at least we are consistent with it, which helps

src/TimeSteppers/clock.jl

src/Units.jl

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

…to ss/rework-fts-1

src/OutputReaders/field_time_series.jl

src/Units.jl

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

src/Units.jl

src/TimeSteppers/clock.jl

glwagner · 2023-09-02T02:41:10Z

We will also have to figure out (perhaps in a future PR) how to insert boundary conditions to Time, eg Cyclic (so forcing / specification repeats) or Continued

glwagner · 2023-09-02T02:48:23Z

src/OutputReaders/field_time_series.jl

+
+# Linear time interpolation
+function Base.getindex(fts::FieldTimeSeries, i::Int, j::Int, k::Int, time_index::Time)
+    Ntimes = length(fts.times)


length is slow for big vectors right? Should we allow this to be passed into Time?

then we can implement something like

how_many_times(time::Time{<:Int}, times) = time.Ntimes how_many_times(::Time, times) = length(times)

or mahybe you have a better name

glwagner · 2023-09-02T02:50:01Z

src/OutputReaders/field_time_series.jl

+    if n₁ == n₂ # no interpolation
+        return getindex(fts, i, j, k, n₁)
+    end
+
+    # fractional index
+    @inbounds n = (n₂ - n₁) / (fts.times[n₂] - fts.times[n₁]) * (time - fts.times[n₁]) + n₁
+    return getindex(fts, i, j, k, n₂) * (n - n₁) + getindex(fts, i, j, k, n₁) * (n₂ - n)
+end


Suggested change

if n₁ == n₂ # no interpolation

return getindex(fts, i, j, k, n₁)

end

# fractional index

@inbounds n = (n₂ - n₁) / (fts.times[n₂] - fts.times[n₁]) * (time - fts.times[n₁]) + n₁

return getindex(fts, i, j, k, n₂) * (n - n₁) + getindex(fts, i, j, k, n₁) * (n₂ - n)

end

# fractional index

@inbounds n = (n₂ - n₁) / (fts.times[n₂] - fts.times[n₁]) * (time - fts.times[n₁]) + n₁

fts_interpolated = getindex(fts, i, j, k, n₂) * (n - n₁) + getindex(fts, i, j, k, n₁) * (n₂ - n)

# Don't interpolate if n = 0.

return ifelse(n₁ == n₂, getindex(fts, i, j, k, n₁), fts_interpolated)

end

I prefer not to do the interpolation if it is not required

Ah, I was actually looking at the other getindex, no in the kernel it's ok

(It's not good practice to resolve other people's comments, unless they are not responding or something --- the reviewer gets priority to decide whether their comment has been resolved or not.)

What kernel is fast? I'm referring to the situation where this is embedded in another kernel. For example, this getindex call can be embedded in the kernels that calculate model tendencies.

On the other hand, a simple if condition should be okay here as there should not be problems of branch divergence (all threads should look at the same time).

Got it. Is that still the case in very complicated kernels, for example kernels that calculate 3D tendencies?

sorry, linear interpolation in GPU kernels is fast. That is if threads access the same time index. There could be a situation in which different threads access different time indexes but it doesn't immediately come to mind.

My question is whether that principle is still useful when the linear interpolation is surrounded by lots of other expensive operations, as can happen when this is used inside a tendency kernel.

Here's another example: FieldTimeSeries is used to store data for an advecting velocity field. This is the application you'd get if you're interested in, for example, advecting biogeochemical tracers with a velocity field derived from the ECCO state estimate. In that case the linear interpolation is inserted into an advection scheme.

* Linear time interpolation * fixed * little bug when n == n1 == n2 * time implementation * Update src/TimeSteppers/clock.jl Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com> * Update src/Units.jl Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com> * fix * Update src/Units.jl Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com> * add triple ticks * import * change --------- Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

simone-silvestri and others added 3 commits August 26, 2023 11:34

Linear time interpolation

82b1abd

fixed

be88fd9

Merge branch 'main' into ss/rework-fts-1

012614d

simone-silvestri added 2 commits August 28, 2023 14:06

little bug when n == n1 == n2

aa39bcb

Merge branch 'ss/rework-fts-1' of github.com:CliMA/Oceananigans.jl in…

82f264f

…to ss/rework-fts-1

time implementation

c535521

glwagner reviewed Aug 29, 2023

View reviewed changes

src/TimeSteppers/clock.jl Outdated Show resolved Hide resolved

glwagner reviewed Aug 29, 2023

View reviewed changes

src/Units.jl Outdated Show resolved Hide resolved

navidcy mentioned this pull request Aug 29, 2023

Better validation for LatitudeLongitudeGrid + add test for sum(Azᶜᶜᵃ) ≈ 4πr² #3240

Merged

simone-silvestri and others added 5 commits August 29, 2023 08:58

Update src/TimeSteppers/clock.jl

25366aa

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

Update src/Units.jl

8ee465d

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

fix

af8f7c6

Merge branch 'ss/rework-fts-1' of github.com:CliMA/Oceananigans.jl in…

30a7df0

…to ss/rework-fts-1

Merge branch 'main' into ss/rework-fts-1

540a21a

glwagner reviewed Aug 29, 2023

View reviewed changes

src/OutputReaders/field_time_series.jl Show resolved Hide resolved

glwagner reviewed Aug 29, 2023

View reviewed changes

src/Units.jl Outdated Show resolved Hide resolved

glwagner approved these changes Aug 29, 2023

View reviewed changes

Update src/Units.jl

98c73e2

Co-authored-by: Gregory L. Wagner <wagner.greg@gmail.com>

glwagner reviewed Aug 29, 2023

View reviewed changes

src/Units.jl Show resolved Hide resolved

add triple ticks

91d86be

glwagner reviewed Aug 30, 2023

View reviewed changes

src/TimeSteppers/clock.jl Show resolved Hide resolved

import

3f14b08

glwagner reviewed Sep 2, 2023

View reviewed changes

change

8a08b17

navidcy added the feature 🌟 Something new and shiny label Sep 7, 2023

simone-silvestri merged commit 2328d54 into main Sep 18, 2023

simone-silvestri deleted the ss/rework-fts-1 branch September 18, 2023 12:19

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Linear time interpolation in `FieldTimeSeries` #3236

Linear time interpolation in `FieldTimeSeries` #3236

simone-silvestri commented Aug 26, 2023 •

edited

Loading

navidcy commented Aug 28, 2023 •

edited

Loading

simone-silvestri commented Aug 28, 2023 •

edited

Loading

navidcy commented Aug 28, 2023

glwagner commented Aug 28, 2023

navidcy commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023 •

edited

Loading

simone-silvestri commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023

glwagner commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023 •

edited

Loading

glwagner commented Aug 28, 2023

glwagner commented Aug 28, 2023

glwagner commented Sep 2, 2023

glwagner Sep 2, 2023

glwagner Sep 2, 2023 •

edited

Loading

glwagner Sep 2, 2023

simone-silvestri Sep 2, 2023

glwagner Sep 3, 2023

simone-silvestri Sep 3, 2023

glwagner Sep 3, 2023

glwagner Sep 3, 2023

glwagner Sep 3, 2023

simone-silvestri Sep 3, 2023

glwagner Sep 4, 2023

glwagner Sep 4, 2023

Linear time interpolation in FieldTimeSeries #3236

Linear time interpolation in FieldTimeSeries #3236

Conversation

simone-silvestri commented Aug 26, 2023 • edited Loading

navidcy commented Aug 28, 2023 • edited Loading

simone-silvestri commented Aug 28, 2023 • edited Loading

navidcy commented Aug 28, 2023

glwagner commented Aug 28, 2023

navidcy commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023 • edited Loading

simone-silvestri commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023

glwagner commented Aug 28, 2023

glwagner commented Aug 28, 2023

simone-silvestri commented Aug 28, 2023 • edited Loading

glwagner commented Aug 28, 2023

glwagner commented Aug 28, 2023

glwagner commented Sep 2, 2023

Choose a reason for hiding this comment

glwagner Sep 2, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Linear time interpolation in `FieldTimeSeries` #3236

Linear time interpolation in `FieldTimeSeries` #3236

simone-silvestri commented Aug 26, 2023 •

edited

Loading

navidcy commented Aug 28, 2023 •

edited

Loading

simone-silvestri commented Aug 28, 2023 •

edited

Loading

simone-silvestri commented Aug 28, 2023 •

edited

Loading

simone-silvestri commented Aug 28, 2023 •

edited

Loading

glwagner Sep 2, 2023 •

edited

Loading