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

Daily interpolation from monthly data (v2) #105

Merged
merged 1 commit into from
Jul 29, 2022
Merged

Daily interpolation from monthly data (v2) #105

merged 1 commit into from
Jul 29, 2022

Conversation

LenkaNovak
Copy link
Collaborator

@LenkaNovak LenkaNovak commented Jul 28, 2022

PULL REQUEST

Purpose and Content

Adds capability to interpolate linearly between two monthly Fields. The demo uses prerequisite PRs outlined in #88 to obtain monthly data.

Benefits and Risks

  • benefit: last piece for enabling specifying temporarily varying boundary conditions from a file.
  • risk: performance slow-down - check are in place (see below)

Linked Issues

Dependent PRs

Performance of the final solution

PR Checklist

  • This PR has a corresponding issue OR is linked to an SDI.
  • I have followed CliMA's codebase contribution and style guidelines OR N/A.
  • I have followed CliMA's documentation policy.
  • I have checked all issues and PRs and I certify that this PR does not duplicate an open PR.
  • I linted my code on my local machine prior to submission OR N/A.
  • Unit tests are included OR N/A.
  • Code used in an integration test OR N/A.
  • All tests ran successfully on my local machine OR N/A.
  • All classes, modules, and function contain docstrings OR N/A.
  • Documentation has been added/updated OR N/A.

@LenkaNovak LenkaNovak marked this pull request as ready for review July 28, 2022 13:31
elseif (midmonth_idx == midmonth_idx0) && (Dates.days(date - all_dates[midmonth_idx]) < 0) # for init
midmonth_idx = bcf_info.segment_idx[1] -= Int(1)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like segment_idx and segment_idx0 are vectors, but the first element is always used. Can these be ints instead of Vectors?

The description of segment_idx0 also makes it sound like it should be a single element.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I wanted to make it a single element, but it needs to be mutable (at least segment_idx needs to be updated every month). Is there a better way of doing this? 🤔 I've heard mutable structs are not great for performance...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good question! Maybe @jb-mackay knows?

I think maybe the answer is that if it is mutating, it should be not stored in in the struct, but instead passed in as an argument to functions that need it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm, not sure exactly how you mean in this context. Would you mind giving an example?

Copy link
Contributor

@kmdeck kmdeck Aug 4, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hey @LenkaNovak ! Happy to chat about this more. This is what I was thinking...

instead of this:

function foo(bc_info::BCDataInfo)
    # do something with bc_info.segment_idx[1]
end

# in main driver couper loop:
 bc_info.segment_idx[1] .= [10] # or however it is updated
foo(bc_info)

we could instead remove segment_idx from the struct, and do something like this

function foo(bc_info::BCDataInfo, segment_idx::Int)
    # does something with segment_idx
end

# in main driver
segment_idx = segment_idx0 # initialize as appropriate
# in coupler loop
segment_idx = 10 # update to next value
foo(bc_info, segment_idx)

Copy link
Contributor

@kmdeck kmdeck Aug 4, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The thought was that if we have an integer valued object that is changing value every step, why store it in the struct? why not update it every step of the main coupler loop?

"""
function interpolate_midmonth_to_daily(date, bcf_info)
@show bcf_info.segment_length[1]
if bcf_info.interpolate_daily[1] && bcf_info.segment_length[1] > FT(0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isnt interpolate_daily a Boolean? why the indexing?

(another comment - in the struct defining BCInfo, we could make interpolate_daily::Bool right away, instead of making it interpolate_daily::I with I a parametric type.)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, this one we can change. Originally I had it set up to be mutable (so that after init it switched to true), but in this version it can be static.

all_dates = bcf_info.all_dates
monthly_fields = bcf_info.monthly_fields

return intepol.(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

intepol -> interpol

"""
intepol(f_t1, f_t2, t, Δt, FT)

Performs linear interpolation within a segment `Δt` = (t1 - t2), of fields `f_t1` and `f_t2`, with `t`(`f_t1`) = t1.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

with f(t_1) = f_1, right? Not t(f_1) = t_1.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe also clearer to have t_2 > t_1, and dt = t2-t1?


Performs linear interpolation within a segment `Δt` = (t1 - t2), of fields `f_t1` and `f_t2`, with `t`(`f_t1`) = t1.
"""
function intepol(f1, f2, t, Δt, FT)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can instead do "interpol(f1::FT, f2::FT, t::FT, dt::FT) where {FT}" and not pass in the FT as well.

"""
function intepol(f1, f2, t, Δt, FT)
interp_fraction = FT(t) / Δt
@assert abs(interp_fraction) <= FT(1) "time interpolation weights must be <= 1, but `interp_fraction` = $interp_fraction"
Copy link
Contributor

@kmdeck kmdeck Jul 28, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dont think this formula is correct.
For example:
If t = t1, we get f1t1/Delta +f2(Delta-t1)/Delta =
f1*t1/Delta - f2 *t2/Delta. It should be f1.

I think it should be f = (f2-f1)/(t2-t1)(t-t1) + f1. Then f(t1) = f1, f(t2) = f2, and f((t1+t2)/2) = (f1+f2)/2.

(Also, it should not matter if t is negative, the formula should work still)

Copy link
Collaborator Author

@LenkaNovak LenkaNovak Jul 28, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ahh, the notation is a little confusing. 🤯 t here actually does refer to difference between the current date and t1, which I believe makes it equivalent to your suggestion. I failed at trying to be less verbose, haha! I will call it \Delta_tt1 and \Delta_t1t2. This is why reviews are super helpful! :) As for the @assert, this is just catching cases where the dates are not passed in properly, e.g., when the current date is outside of range [t1,t2].

segment accuracy fixed

test update

rev fixes
@LenkaNovak
Copy link
Collaborator Author

LenkaNovak commented Jul 29, 2022

All issues above have been addressed, except whether to use functions instead of the mutable vectors in static structs. At the last sync meeting it was decided that since this is not a breaking issue, we will merge this and address this as part of a future PR. This was recorded as a TODO in issue #92

@LenkaNovak
Copy link
Collaborator Author

bors r+

@bors
Copy link
Contributor

bors bot commented Jul 29, 2022

@bors bors bot merged commit ec8de33 into main Jul 29, 2022
@bors bors bot deleted the daily-interpol_v2 branch July 29, 2022 21:21
@LenkaNovak LenkaNovak mentioned this pull request Sep 7, 2022
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants