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.

Sign up for GitHub

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

Jump to bottom

Fix initial file read #579

Merged
merged 1 commit into from
Jan 31, 2024
Merged

Fix initial file read #579

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
100 changes: 63 additions & 37 deletions src/BCReader.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ and returns the info packaged in a single struct.
- `boundary_space`: [Spaces.AbstractSpace] the space to which we are mapping.
- `comms_ctx`: [ClimaComms.AbstractCommsContext] context used for this operation.
- `interpolate_daily`: [Bool] switch to trigger daily interpolation.
- `segment_idx0`: [Vector{Int}] index of the file data that is closest to date0.
- `segment_idx0`: [Vector{Int}] reference date which, after initialization, refers to the the first file date index used minus 1 (segment_idx[1] - 1)
- `scaling function`: [Function] scales, offsets or transforms `varname`.
- `land_fraction`: [Fields.field] fraction with 1 = land, 0 = ocean / sea-ice.
- `date0`: [Dates.DateTime] start date of the file data.
Expand Down Expand Up @@ -182,22 +182,54 @@ The times for which data is extracted depends on the specifications in the
function update_midmonth_data!(date, bcf_info::BCFileInfo{FT}) where {FT}
# monthly count
(; bcfile_dir, comms_ctx, hd_outfile_root, varname, all_dates, scaling_function) = bcf_info
midmonth_idx = bcf_info.segment_idx[1]
midmonth_idx0 = bcf_info.segment_idx0[1]

if (midmonth_idx == midmonth_idx0) && (Dates.days(date - all_dates[midmonth_idx]) < 0) # for init
midmonth_idx = bcf_info.segment_idx[1] -= Int(1)
midmonth_idx = midmonth_idx < Int(1) ? midmonth_idx + Int(1) : midmonth_idx
@warn "this time period is before BC data - using file from $(all_dates[midmonth_idx0])"
bcf_info.monthly_fields[1] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx0)], varname, comms_ctx),
bcf_info,
)
bcf_info.monthly_fields[2] .= deepcopy(bcf_info.monthly_fields[1])
bcf_info.segment_length .= Int(0)

elseif Dates.days(date - all_dates[end - 1]) > 0 # for fini
@warn "this time period is after BC data - using file from $(all_dates[end - 1])"
segment_idx = bcf_info.segment_idx[1] # index of the current date in the file. [segment_idx, segment_idx+1] indexes the current segment between which we interpolate
segment_idx0 = bcf_info.segment_idx0[1] # reference index (first segment_idx - 1)

# upon initialization (segment_idx == segment_idx0) with model date before the final file date
if (segment_idx == segment_idx0) && !((date - all_dates[end]).value >= 0)
# case 1: model date is before the first segment read from file
if (date - all_dates[segment_idx0]).value < 0
@warn "This time period is before BC data - using file from $(all_dates[segment_idx0])"
bcf_info.monthly_fields[1] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(segment_idx0)], varname, comms_ctx),
bcf_info,
)
bcf_info.monthly_fields[2] .= deepcopy(bcf_info.monthly_fields[1])
bcf_info.segment_length .= Int(0)
bcf_info.segment_idx[1] -= Int(1)
bcf_info.segment_idx0[1] -= Int(2)

# case 2: model date is after the first segment read from file
elseif (date - all_dates[Int(segment_idx0) + 1]).value >= 0
nearest_idx = argmin(
abs.(
parse(FT, TimeManager.datetime_to_strdate(date)) .-
parse.(FT, TimeManager.datetime_to_strdate.(all_dates[:]))
),
)
@warn "Initializing with `segment_idx = $nearest_idx"
bcf_info.segment_idx[1] = nearest_idx
bcf_info.segment_idx0[1] = nearest_idx
update_midmonth_data!(date, bcf_info)

# case 3: model date is within the first segment read from file
LenkaNovak marked this conversation as resolved.
Show resolved Hide resolved
elseif (date - all_dates[segment_idx0]).value >= 0
@warn "On $date updating file data reads: file dates = [ $(all_dates[segment_idx]) , $(all_dates[segment_idx+1]) ]"
LenkaNovak marked this conversation as resolved.
Show resolved Hide resolved
bcf_info.segment_length .= (all_dates[segment_idx + 1] - all_dates[segment_idx]).value
bcf_info.monthly_fields[1] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[segment_idx], varname, comms_ctx),
bcf_info,
)
bcf_info.monthly_fields[2] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[segment_idx + 1], varname, comms_ctx),
bcf_info,
)
bcf_info.segment_idx0[1] -= Int(1)
LenkaNovak marked this conversation as resolved.
Show resolved Hide resolved
end

# case 4: date is at or after the last date in file
elseif (date - all_dates[end]).value >= 0
@warn "This time period is after BC data - using file from $(all_dates[end])"
bcf_info.monthly_fields[1] .= scaling_function(
Regridder.read_from_hdf5(
bcfile_dir,
Expand All @@ -211,34 +243,28 @@ function update_midmonth_data!(date, bcf_info::BCFileInfo{FT}) where {FT}
bcf_info.monthly_fields[2] .= deepcopy(bcf_info.monthly_fields[1])
bcf_info.segment_length .= Int(0)

# throw error when there are closer initial indices for the bc file data that matches this date0
elseif Dates.days(date - all_dates[Int(midmonth_idx + 1)]) > 2
nearest_idx = argmin(
abs.(
parse(FT, TimeManager.datetime_to_strdate(date)) .-
parse.(FT, TimeManager.datetime_to_strdate.(all_dates[:]))
),
)
# TODO test this
bcf_info.segment_idx[1] = midmonth_idx = midmonth_idx0 = nearest_idx
@warn "init data does not correspond to start date. Initializing with `SIC_info.segment_idx[1] = midmonth_idx = midmonth_idx0 = $nearest_idx` for this start date"

# date crosses to the next month
elseif Dates.days(date - all_dates[Int(midmonth_idx)]) > 0
midmonth_idx = bcf_info.segment_idx[1] += Int(1)
@warn "On $date updating monthly data files: mid-month dates = [ $(all_dates[Int(midmonth_idx)]) , $(all_dates[Int(midmonth_idx+1)]) ]"
bcf_info.segment_length .= (all_dates[Int(midmonth_idx + 1)] - all_dates[Int(midmonth_idx)]).value
# case 5: model date crosses to the next segment
elseif (date - all_dates[Int(segment_idx) + 1]).value >= 0
segment_idx = bcf_info.segment_idx[1] += Int(1)

bcf_info.segment_length .= (all_dates[Int(segment_idx + 1)] - all_dates[Int(segment_idx)]).value

bcf_info.monthly_fields[1] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx)], varname, comms_ctx),
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(segment_idx)], varname, comms_ctx),
bcf_info,
)
bcf_info.monthly_fields[2] .= scaling_function(
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(midmonth_idx + 1)], varname, comms_ctx),
Regridder.read_from_hdf5(bcfile_dir, hd_outfile_root, all_dates[Int(segment_idx + 1)], varname, comms_ctx),
bcf_info,
)

# case 6: undefined condition
else
throw(ErrorException("Check boundary file specification"))
throw(
ErrorException(
"Check boundary file specification: segment: $(all_dates[segment_idx]) - $(all_dates[segment_idx+1]), date: $date",
),
)
end
end

Expand Down
131 changes: 92 additions & 39 deletions test/bcreader_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -227,60 +227,113 @@ for FT in (Float32, Float64)
bcf_info.monthly_fields[1] .= current_fields[1]
bcf_info.monthly_fields[2] .= current_fields[2]
bcf_info.segment_length[1] = Int(1)
bcf_info.segment_idx[1] = bcf_info.segment_idx0[1]
end

hd_outfile_root = varname * "_cgll"

# case 1: segment_idx == segment_idx0, date < all_dates[segment_idx]
# case 1: date < all_dates[segment_idx] (init)
bcf_info.segment_idx[1] = bcf_info.segment_idx0[1]
date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1))
BCReader.update_midmonth_data!(date, bcf_info)

@test bcf_info.monthly_fields[1] == bcf_info.scaling_function(
Regridder.read_from_hdf5(
regrid_dir,
hd_outfile_root,
bcf_info.all_dates[Int(bcf_info.segment_idx0[1])],
varname,
comms_ctx,
),
bcf_info,
)
# unmodified field
@test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1]
# zero segment length
@test bcf_info.segment_length[1] == Int(0)
# segment index is reset
@test bcf_info.segment_idx0[1] == bcf_info.segment_idx[1] - 1

# cases 2 and 3
extra_days = [Dates.Day(0), Dates.Day(3)]
for extra in extra_days
# case 3: (date - all_dates[Int(segment_idx0)]) >= 0 (init)
reset_bcf_info(bcf_info)
date = DateTime(bcf_info.all_dates[bcf_info.segment_idx0[1]]) + extra
BCReader.update_midmonth_data!(date, bcf_info)

end_field_c2 = deepcopy(bcf_info.monthly_fields[2])
segment_length_c2 = deepcopy(bcf_info.segment_length[1])
current_index_c2 = deepcopy(bcf_info.segment_idx[1])

# modified field
@test end_field_c2 !== bcf_info.monthly_fields[1]
# updated segment length
@test segment_length_c2[1] !== Int(0)
# updated reference segment index
@test current_index_c2 == bcf_info.segment_idx0[1] + 1

# case 2: (date - all_dates[Int(segment_idx0) + 1]) >= 0 (init)
# do not reset segment_idx0. It's current value ensures that we get the same result as case 3
reset_bcf_info(bcf_info)

date = DateTime(bcf_info.all_dates[bcf_info.segment_idx0[1] + 1]) + extra
BCReader.update_midmonth_data!(date, bcf_info)

nearest_idx = argmin(
abs.(
parse(FT, TimeManager.datetime_to_strdate(date)) .-
parse.(FT, TimeManager.datetime_to_strdate.(bcf_info.all_dates[:]))
),
)

@test bcf_info.segment_idx[1] == bcf_info.segment_idx0[1] + 1 == nearest_idx

# compare to case 3 (expecting the same result - this defaults to it):
@test bcf_info.monthly_fields[1] == bcf_info.scaling_function(
Regridder.read_from_hdf5(
regrid_dir,
hd_outfile_root,
bcf_info.all_dates[Int(bcf_info.segment_idx[1])],
varname,
comms_ctx,
),
bcf_info,
)

# check case 2 defaults to case 3
@test end_field_c2 !== bcf_info.monthly_fields[1]
@test segment_length_c2[1] !== Int(0)
@test current_index_c2 == bcf_info.segment_idx0[1] + 1

# case 2: date > all_dates[end - 1]
reset_bcf_info(bcf_info)
date = DateTime(bcf_info.all_dates[end - 1] + Dates.Day(1))
BCReader.update_midmonth_data!(date, bcf_info)
end

@test bcf_info.monthly_fields[1] == bcf_info.scaling_function(
Regridder.read_from_hdf5(
regrid_dir,
hd_outfile_root,
bcf_info.all_dates[Int(length(bcf_info.all_dates))],
varname,
comms_ctx,
),
bcf_info,
)
@test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1]
@test bcf_info.segment_length[1] == Int(0)
# case 4: date > all_dates[end]
for extra in extra_days
bcf_info.segment_idx0[1] = length(bcf_info.all_dates)
reset_bcf_info(bcf_info)

date = DateTime(bcf_info.all_dates[bcf_info.segment_idx0[1]]) + extra
BCReader.update_midmonth_data!(date, bcf_info)

@test bcf_info.monthly_fields[1] == bcf_info.scaling_function(
Regridder.read_from_hdf5(
regrid_dir,
hd_outfile_root,
bcf_info.all_dates[Int(length(bcf_info.all_dates))],
varname,
comms_ctx,
),
bcf_info,
)
@test bcf_info.monthly_fields[2] == bcf_info.monthly_fields[1]
@test bcf_info.segment_length[1] == Int(0)
end

# case 3: date - all_dates[segment_idx + 1] > 2
reset_bcf_info(bcf_info)
date = DateTime(bcf_info.all_dates[bcf_info.segment_idx[1] + 1] + Dates.Day(3))
BCReader.update_midmonth_data!(date, bcf_info)
# case 5: Dates.days(date - all_dates[segment_idx]) >= 0

nearest_idx = argmin(
abs.(
parse(FT, TimeManager.datetime_to_strdate(date)) .-
parse.(FT, TimeManager.datetime_to_strdate.(bcf_info.all_dates[:]))
),
)
@test bcf_info.segment_idx[1] == bcf_info.segment_idx0[1] == nearest_idx
extra = Dates.Day(3)
for extra in extra_days
bcf_info.segment_idx0[1] = 2
reset_bcf_info(bcf_info)

date = DateTime(bcf_info.all_dates[bcf_info.segment_idx0[1]] + extra)
BCReader.update_midmonth_data!(date, bcf_info)

@test bcf_info.segment_idx[1] == bcf_info.segment_idx0[1] + 1
end

# case 4: everything else
# case 6: everything else
reset_bcf_info(bcf_info)
bcf_info.segment_idx[1] = bcf_info.segment_idx0[1] + Int(1)
date = bcf_info.all_dates[bcf_info.segment_idx[1]] - Dates.Day(1)
Expand Down
Loading
Loading