-
Notifications
You must be signed in to change notification settings - Fork 402
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
improve floating point comparisons in vic_init.c #701
Conversation
@@ -986,8 +986,7 @@ vic_init(void) | |||
} | |||
sum += soil_con[i].AreaFract[j]; | |||
} | |||
// TBD: Need better check for equal to 1. | |||
if (sum != 1.) { | |||
if (!assert_close_double(sum, 1.0, 0., 0.001)) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It makes sense to not throw out a warning for very small sum errors. Is 0.001 enough a cutoff for not adjusting AreaFract though? (Similar for a couple of other changes below)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's have @tbohn and @bartnijssen weigh in here.
Options are:
- Always divide by the sum and only print a warning when the difference is larger than a threshold (0.001).
- Only divide by the sum when the difference is larger than a very small threshold (1e-20), print a warning in this case.
- Raise an error in cases where the difference is greater than a very small threshold (1e-20) instead of doing work for the user.
Thoughts?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess option 3 allows for the simplest, easiest-to-maintain code. But it does place the burden on the user - is it possible to have a case where a user could have real trouble getting fractions to sum to within 1e-20 of 1.0 in their netcdf file? Or perhaps a better question is: would it be common enough for this to happen that VIC should be more accommodating?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Generally I prefer to have the work done by the user since that is the most explicit. They should only need to do this once and it seems to me that that would be an easy calculation in some version of nco
, cdo
or xarray
. 0.001
would seem too large to me to take care of behind the scenes. The one thing is that our documentation needs to be very clear on what we expect from the user, so if we print an error message, then our documentation needs to be explicit on how to fix it or prevent it.
That said, I can live with 2.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks @bartnijssen and @tbohn. I will update this PR to include a smaller threshold (probably 1e-10) and add clear warnings/documentation to that effect.
@UW-Hydro/vic-developers - this is ready for another review. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll wait with further review till it's resolved that this is what we want
@@ -42,6 +42,7 @@ calc_root_fractions(veg_con_struct *veg_con, | |||
size_t layer; | |||
size_t zone; | |||
size_t i; | |||
size_t max_itter; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
should read max_iter
. Actually it would make more sense if this were simply called iter
or n_iter
, since it is not the maximum number of iterations but just a counter
while (zone < options.ROOT_ZONES) { | ||
max_itter++; | ||
if (max_itter > 9999) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Replace 9999
by some #define
or constant, e.g. MAX_ITER
. Any reason why this is set to 9999
that seems a large number of iterations.
@@ -62,7 +63,16 @@ calc_root_fractions(veg_con_struct *veg_con, | |||
Zsum = 0; | |||
zone = 0; | |||
|
|||
max_itter = -1; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why not start at 0
and update the iterator at the end of the iteration rather than the beginning?
@@ -124,6 +134,7 @@ calc_root_fractions(veg_con_struct *veg_con, | |||
} | |||
} | |||
else if (Zsum + Zstep > Lsum) { | |||
zone++; // should this be here? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No idea - shouldn't that be up to the person changing the code? At the very least the comment should be resolved (and removed) before the PR is accepted.
…t_library_shared_image.c
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Same comments as earlier. If you don't change them - just reply to the comment so I know you looked at it.
@@ -124,6 +134,7 @@ calc_root_fractions(veg_con_struct *veg_con, | |||
} | |||
} | |||
else if (Zsum + Zstep > Lsum) { | |||
zone++; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@tbohn - can you look at the changes in this file, specifically this line?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good. I'm assuming you've tested this to ensure that the root fractions come out the same as before?
I wonder if there's a simpler algorithm for this function. Something to look into some day.
@bartnijssen - I think we're all good here. Want to take one more look? |
@jhamman - review passed - just merging develop back into it since it said it was out-of-date. If the check passes I'll go ahead and merge. |
A few cleanup items.