Skip to content
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

Changes the internal behavior of TimeVaryingInputs23D to use dates rather than seconds since start date #147

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/src/datahandling.md
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,8 @@ ClimaUtilities.DataHandling.available_times
ClimaUtilities.DataHandling.available_dates
ClimaUtilities.DataHandling.previous_time
ClimaUtilities.DataHandling.next_time
ClimaUtilities.DataHandling.previous_date
ClimaUtilities.DataHandling.next_date
ClimaUtilities.DataHandling.regridded_snapshot
ClimaUtilities.DataHandling.regridded_snapshot!
ClimaUtilities.DataHandling.dt
Expand Down
40 changes: 40 additions & 0 deletions ext/DataHandlingExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -511,6 +511,46 @@ function DataHandling.next_time(data_handler::DataHandler, date::Dates.DateTime)
return data_handler.available_times[index]
end

"""
DataHandling.previous_date(data_handler::DataHandler, time::Dates.TimeType)

Return the date of the snapshot before the given `date`.
If `date` is one of the snapshots, return itself.

If `date` is not in the `data_handler`, return an error.
"""
function DataHandling.previous_date(
data_handler::DataHandler,
date::Dates.TimeType,
)
if date in data_handler.available_dates
index = searchsortedfirst(data_handler.available_dates, date)
else
index = searchsortedfirst(data_handler.available_dates, date) - 1
end
index < 1 && error("Date $date is before available dates")
return data_handler.available_dates[index]
end

"""
DataHandling.next_date(data_handler::DataHandler, time::Dates.TimeType)

Return the date of the snapshot after the given `time`.
If `date` is one of the snapshots, return itself.

If `date` is not in the `data_handler`, return an error.
"""
function DataHandling.next_date(data_handler::DataHandler, date::Dates.TimeType)
if date in data_handler.available_dates
index = searchsortedfirst(data_handler.available_dates, date) + 1
else
index = searchsortedfirst(data_handler.available_dates, date)
end
index > length(data_handler.available_dates) &&
error("Date $date is after available dates")
return data_handler.available_dates[index]
end

"""
regridded_snapshot(data_handler::DataHandler, date::Dates.DateTime)
regridded_snapshot(data_handler::DataHandler, time::AbstractFloat)
Expand Down
90 changes: 58 additions & 32 deletions ext/TimeVaryingInputsExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,12 +33,13 @@ import ClimaUtilities.TimeVaryingInputs: extrapolation_bc
import ClimaUtilities.DataHandling
import ClimaUtilities.DataHandling:
DataHandler,
previous_time,
next_time,
regridded_snapshot!,
available_times,
available_dates,
date_to_time
date_to_time,
time_to_date,
previous_date,
next_date

# Ideally, we should be able to split off the analytic part in a different
# extension, but precompilation stops working when we do so
Expand Down Expand Up @@ -110,9 +111,16 @@ end
Check if the given `time` is in the range of definition for `itp`.
"""
function Base.in(time, itp::InterpolatingTimeVaryingInput23D)
return itp.range[1] <= time <= itp.range[2]
return itp.data_handler.available_dates[begin] <=
time <=
itp.data_handler.available_dates[end]
end

function Base.in(time::Number, itp::InterpolatingTimeVaryingInput23D)
return Base.in(time_to_date(itp.data_handler, time), itp)
end


function TimeVaryingInputs.TimeVaryingInput(
data_handler;
method = LinearInterpolation(),
Expand Down Expand Up @@ -198,19 +206,37 @@ function TimeVaryingInputs.evaluate!(
time in itp || error("TimeVaryingInput does not cover time $time")
end
if extrapolation_bc(itp.method) isa Flat
t_init, t_end = itp.range
if time >= t_end
regridded_snapshot!(dest, itp.data_handler, t_end)
date_init, date_end = itp.data_handler.available_dates[begin],
itp.data_handler.available_dates[end]
if time >= date_end
regridded_snapshot!(dest, itp.data_handler, date_end)
else
time <= t_init
regridded_snapshot!(dest, itp.data_handler, t_init)
time <= date_init
regridded_snapshot!(dest, itp.data_handler, date_init)
end
else
TimeVaryingInputs.evaluate!(dest, itp, time, itp.method)
end
return nothing
end

function TimeVaryingInputs.evaluate!(
dest,
itp::InterpolatingTimeVaryingInput23D,
time::Number,
args...;
kwargs...,
)
TimeVaryingInputs.evaluate!(
dest,
itp,
Dates.Millisecond(round(1_000 * time)) + itp.data_handler.start_date,
args...,
kwargs...,
)
return nothing
end

function _time_range_dt_dt_e(itp::InterpolatingTimeVaryingInput23D)
return _time_range_dt_dt_e(itp, extrapolation_bc(itp.method))
end
Expand All @@ -220,7 +246,10 @@ function _time_range_dt_dt_e(
extrapolation_bc::PeriodicCalendar{Nothing},
)
dt = DataHandling.dt(itp.data_handler)
return itp.range[begin], itp.range[end], dt, dt / 2
return itp.data_handler.available_dates[begin],
itp.data_handler.available_dates[end],
Dates.Millisecond(round(1_000 * dt)),
Dates.Millisecond(round((1_000 * dt) / 2))
end

function _time_range_dt_dt_e(
Expand All @@ -240,11 +269,9 @@ function _time_range_dt_dt_e(
date_to_time(itp.data_handler, date_end)
# We have to add 1 Second because endofperiod(date_end, period) returns the very last
# second before the next period
dt_e =
(endofperiod(date_end, period) + Dates.Second(1) - date_end) |>
period_to_seconds_float
dt = (date_init + period - date_end) |> period_to_seconds_float
return t_init, t_end, dt, dt_e
dt_e = (endofperiod(date_end, period) + Dates.Second(1) - date_end)
dt = (date_init + period - date_end)
return date_init, date_end, dt, dt_e
end

"""
Expand Down Expand Up @@ -320,15 +347,14 @@ function TimeVaryingInputs.evaluate!(
end
return nothing
end

t0, t1 =
previous_time(itp.data_handler, time), next_time(itp.data_handler, time)
date0, date1 =
previous_date(itp.data_handler, time), next_date(itp.data_handler, time)

# The closest regridded_snapshot could be either the previous or the next one
if (time - t0) <= (t1 - time)
regridded_snapshot!(dest, itp.data_handler, t0)
if (time - date0) <= (date1 - time)
regridded_snapshot!(dest, itp.data_handler, date0)
else
regridded_snapshot!(dest, itp.data_handler, t1)
regridded_snapshot!(dest, itp.data_handler, date1)
end
return nothing
end
Expand All @@ -351,7 +377,7 @@ function TimeVaryingInputs.evaluate!(
field_t0, field_t1 = itp.preallocated_regridded_fields[begin:(begin + 1)]

if extrapolation_bc(itp.method) isa PeriodicCalendar
time, t_init, t_end, dt, _ =
time, date_init, date_end, dt, _ =
_interpolation_times_periodic_calendar(time, itp)

# We have to handle separately the edge case where the desired time is past t_end.
Expand All @@ -361,26 +387,26 @@ function TimeVaryingInputs.evaluate!(

# TODO: It would be nice to handle this edge case directly instead of copying the
# code
if time > t_end
regridded_snapshot!(field_t0, itp.data_handler, t_end)
regridded_snapshot!(field_t1, itp.data_handler, t_init)
coeff = (time - t_end) / dt
if time > date_end
regridded_snapshot!(field_t0, itp.data_handler, date_end)
regridded_snapshot!(field_t1, itp.data_handler, date_init)
coeff = (time - date_end) / dt
dest .= (1 - coeff) .* field_t0 .+ coeff .* field_t1
return nothing
end
end

# We have to consider the edge case where time is precisely the last available_time.
# This is relevant also because it can be triggered by LinearPeriodFilling
if time in DataHandling.available_times(itp.data_handler)
if time in DataHandling.available_dates(itp.data_handler)
regridded_snapshot!(dest, itp.data_handler, time)
else
t0, t1 = previous_time(itp.data_handler, time),
next_time(itp.data_handler, time)
coeff = (time - t0) / (t1 - t0)
date0, date1 = previous_date(itp.data_handler, time),
next_date(itp.data_handler, time)
coeff = (time - date0) / (date1 - date0)

regridded_snapshot!(field_t0, itp.data_handler, t0)
regridded_snapshot!(field_t1, itp.data_handler, t1)
regridded_snapshot!(field_t0, itp.data_handler, date0)
regridded_snapshot!(field_t1, itp.data_handler, date1)

dest .= (1 - coeff) .* field_t0 .+ coeff .* field_t1
end
Expand Down
35 changes: 24 additions & 11 deletions ext/time_varying_inputs_linearperiodfilling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ function TimeVaryingInputs.evaluate!(
available_periods = unique_periods(available_dates, period)

# E.g., Date(1987, 2, 1)
target_date = DataHandling.time_to_date(itp.data_handler, time)
target_date = time

# If the target date falls within an available period, we just interpolate it with
# LinearInterpolation()
Expand Down Expand Up @@ -274,8 +274,6 @@ function TimeVaryingInputs.evaluate!(
if in_interpolable_region
date_pre = _move_date_to_period(target_date, period_left, period)
date_post = _move_date_to_period(target_date, period_right, period)
time_pre = date_to_time(itp.data_handler, date_pre)
time_post = date_to_time(itp.data_handler, date_post)

# Linear interpolation: y = y0 * (1 - coeff) + coeff * y1
#
Expand All @@ -288,12 +286,12 @@ function TimeVaryingInputs.evaluate!(
period_offset = period_right - period_left
coeff = offset_periods_left / period_offset

TimeVaryingInputs.evaluate!(dest, itp, time_pre, LinearInterpolation())
TimeVaryingInputs.evaluate!(dest, itp, date_pre, LinearInterpolation())
dest .*= (1 - coeff)
TimeVaryingInputs.evaluate!(
tmp_field1,
itp,
time_post,
date_post,
LinearInterpolation(),
)
tmp_field1 .*= coeff
Expand Down Expand Up @@ -359,17 +357,32 @@ function TimeVaryingInputs.evaluate!(
error("We should not be here!")
end

time_pre = date_to_time(itp.data_handler, date_pre)
time_post = date_to_time(itp.data_handler, date_post)

# y = y0 * (1 - coeff) + coff * y1
TimeVaryingInputs.evaluate!(dest, itp, time_pre, method)
coeff = (time - time_pre) / (time_post - time_pre)
TimeVaryingInputs.evaluate!(dest, itp, date_pre, method)
coeff = (time - date_pre) / (date_post - date_pre)
dest .*= (1 - coeff)
TimeVaryingInputs.evaluate!(tmp_field2, itp, time_post, method)
TimeVaryingInputs.evaluate!(tmp_field2, itp, date_post, method)
tmp_field2 .*= coeff
dest .+= tmp_field2
return nothing
end
return nothing
end

function TimeVaryingInputs.evaluate!(
dest,
itp::InterpolatingTimeVaryingInput23D,
time::Number,
method::LinearPeriodFillingInterpolation,
args...;
kwargs...,
)
TimeVaryingInputs.evaluate!(
dest,
itp,
Dates.Millisecond(round(1_000 * time)) + itp.data_handler.start_date,
args...,
kwargs...,
)
return nothing
end
8 changes: 8 additions & 0 deletions src/DataHandling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,10 @@ function previous_time end

function next_time end

function previous_date end

function next_date end

function regridded_snapshot end

function regridded_snapshot! end
Expand All @@ -68,6 +72,8 @@ extension_fns = [
:available_dates,
:previous_time,
:next_time,
:previous_date,
:next_date,
:regridded_snapshot,
:regridded_snapshot!,
:dt,
Expand All @@ -80,6 +86,8 @@ extension_fns = [
:available_dates,
:previous_time,
:next_time,
:previous_date,
:next_date,
:regridded_snapshot,
:regridded_snapshot!,
:dt,
Expand Down
50 changes: 50 additions & 0 deletions test/data_handling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -361,6 +361,31 @@ end
available_dates[10],
) == available_times[10]

# Previous date
@test_throws ErrorException DataHandling.previous_date(
data_handler,
available_dates[1] - Second(1),
)

@test DataHandling.previous_date(
data_handler,
available_dates[10] + Second(1),
) == available_dates[10]
@test DataHandling.previous_date(
data_handler,
available_dates[1] + Second(1),
) == available_dates[1]

# Previous time with time, boundaries (return the node)
@test DataHandling.previous_date(
data_handler,
available_dates[1],
) == available_dates[1]
@test DataHandling.previous_date(
data_handler,
available_dates[end],
) == available_dates[end]

# Next time with time
@test_throws ErrorException DataHandling.next_time(
data_handler,
Expand Down Expand Up @@ -391,6 +416,31 @@ end
available_dates[10],
) == available_times[11]

# Next date
@test_throws ErrorException DataHandling.next_date(
data_handler,
available_dates[end] + Second(1),
)

@test DataHandling.next_date(
data_handler,
available_dates[10] + Second(1),
) == available_dates[11]
@test DataHandling.next_date(
data_handler,
available_dates[1] + Second(1),
) == available_dates[2]

# On node
@test DataHandling.next_date(
data_handler,
available_dates[10],
) == available_dates[11]
@test DataHandling.next_date(
data_handler,
available_dates[10],
) == available_dates[11]

# Asking for a regridded_snapshot without specifying the time
@test_throws ErrorException DataHandling.regridded_snapshot(
data_handler,
Expand Down
Loading
Loading