Comments (66)
Are you gathering data using the same project as for #10?
from crhmcode.
Pretty much, although I used a SWE threshold of 1 mm.
This one calcualted the radiation components
radiation_model.zip
This one had calculaed the net radiation
Evap_model.zip
from crhmcode.
What HRU's are you seeing this difference in?
from crhmcode.
I was looking at HRU 11
from crhmcode.
I am not seeing differences in net, netD, or RnD that are outside of our tolerance of 5% or 0.0001.
I checked hru's 2, 3, 11 for the first two years of the simulation.
What years are you examining?
from crhmcode.
Sorry, the examples I plotted above are from 1995.
You can see that in many years, we get very good agreement between both models.
However in some years, chrmcode is underestimating runoff, probably because the evaporation
is too large
from crhmcode.
And you believe that the evaporation is too large because our radiation values are too large as well?
from crhmcode.
Probably - haven't checked it yet, but it's the most likely scenario.
I'll do a more complete water balance to be sure.
from crhmcode.
Alright I will investigate the overflow, evap, and rad, variables around 1960-1968 and then see if I can figure out what is causing it.
from crhmcode.
I plotted the daily evap for HRU 11 over the whole run period.
Most of the values are very similar, but there are quite a few outliers
The crhmcode mean daily evap is slightly smaller than the Borland evap for HRU 11
date Borland crhmcode
Min. :1960-01-01 Min. :0.0000 Min. :0.0000
1st Qu.:1971-08-06 1st Qu.:0.0000 1st Qu.:0.0000
Median :1983-03-11 Median :0.7606 Median :0.7227
Mean :1983-03-11 Mean :2.1168 Mean :2.1071
3rd Qu.:1994-10-14 3rd Qu.:4.4962 3rd Qu.:4.4766
Max. :2006-05-20 Max. :8.0753 Max. :8.0753
from crhmcode.
I have found that the albedo is not matching for hru 11 I assume that this will affect the evap?
from crhmcode.
yes it will - this is important.
from crhmcode.
I got Logan to output the values of some of the internal variables from the Borland netall for HRU 11
for May 21-28, 1995.
netall_HRU11_May_21-28_1995_run2.txt
from crhmcode.
I output the daily evaporation, radiation and albedo for days where the Broaldn evaporation was non zero and > 1.1 x the crhmcode evaporation and vice versa.
The biggest thing that stands out is the number of days where one of the models was setting the albedo to 0.11, i.e. to soil rather than snow.
HRU11_daily_evap_netrad_albedo.zip
I wonder if this might be due to a differing precision when evaluating
if(SWE[hh] > 5.0) {
Albedo[hh] = Albedo_snow[hh];
winter[hh] = 1;
}
else{
Albedo[hh] = Albedo_bare[hh];
winter[hh] = 0;
}
in classalbedo
from crhmcode.
That certainly is possible. Anytime we have comparisons with floating points we are going to have issues as comparisons on floating points are never exact.
Especially with moving from single to double precision.
from crhmcode.
I implemented the instrumentation that Logan did and analyzed the same time period.
instrumented_results.zip
Results reveal that the instrumented variables are identical save for precision error on that time period but the radiation values do not match.
from crhmcode.
I ran a test for the year of 1995 to see if that swe comparison was not being evaluated the same. I used the winter varaiable as a proxy to see what branch was evaluated.
For that year I did not see any times where borland and crhmcode differed I will try with a longer timeline on Monday.
from crhmcode.
I think I might have it.
I output meltflag for both models and crhmcode is seting meltflag to 1 in the interval 1995-06-21 00:00 through
1995-06-26 23:00, while the Borland version doesn't. It's being set by this condition which also sets the
albedo:
if(SWE[hh] <= 0.0) {
Albedo[hh] = Albedo_bare[hh];
winter[hh] = 0;
meltflag[hh] = 0;
}
Looking at the SWE values, the crhmcode values are, again, very slightly greater than zero,
while the Borland versions are not.
We need to use a threshold for the SWE test. Why don't we try using the existing SWE threshold?
from crhmcode.
I think we should also check out the ratio of rainfall to precipitation calculated in classobs
hru_rain[hh] = hru_p[hh]*ratio;
hru_snow[hh] = hru_p[hh]*(1.0-ratio);
We are getting some very small values of SWE, which really should be taken care of
If we set the ratio to 1.0 and 0.0 for values >= 0.999 and <= 0.001 (or similar values) then we should avoid the roundoff error SWE values
from crhmcode.
We need to use a threshold for the SWE test. Why don't we try using the existing SWE threshold?
We can certainly do this. My only concern would be that the parameter as we had specified it was a specific threshold for preventing evaporation. If you feel using the same number here makes sense in terms of the physical process we can certainly apply it.
We will have to get Logan to change the newly added parameter to a shared parameter.
However I think your second point is more important. If we can fix the root of the very small amount of snow problem this would possibly be the best in general. I will see if I can come up with reasonable numbers that have the desired effect.
from crhmcode.
Those are very good points. I agree that it's best to fix the root of the problem and to keep the code as similar as possible between the two versions
from crhmcode.
I added some code so I could see what ratio was set to at various points and looked at the first two years of the model.
As you observed there are many points where this is leading to extremely small values for rain or snow. However, simply setting the ratio to 0 or 1 when it approaches that value will not solve our problem. As the data I posted shows there are points where the ratio value is set as 1 but the nature of floating points and the fact that snow is calculated as ppt* (1.0 - ratio) means that we get extremely small amounts of snow.
I will see if I can come up with a more comprehensive fix.
from crhmcode.
I was afraid of that. What would happen if we just set snow = 0.0 or snow = 1.0, instead of setting the ratio?
Would the roundoff error still occur?
from crhmcode.
I think that is a more correct approach and what I was going to explore first.
hru_snow[hh] = 0.0;
hru_rain[hh] = 0.0;
if(hru_p[hh] > 0.0) //rain or snow determined by ice bulb ratio
hru_rain[hh] = hru_p[hh]*ratio;
hru_snow[hh] = hru_p[hh]*(1.0-ratio);
In a related area the above code is the calucualtion we are discussing and I think it might be having other unintended behavior but I would like to run it by you.
As it is written now when hru_p is less than zero we are guaranteed to have a value of 0 for rain however their can still be a non zero value for snow is this intended?
from crhmcode.
No, both rain and snow should never be less than zero. More bad coding (although IIRC the values are checked somewhere else).
from crhmcode.
But more specifically if hru_p is 0.0 then both hru_rain and hru_snow should be 0.0 as well?
Is it also safe to assume that hru_p = hru_rain + hru_snow?
from crhmcode.
Exactly. Conservation of mass requires that rain + snow = p.
Neither can ever be negative.
from crhmcode.
Alright if that is the desired outcome that if statement should have both the setting of rain and snow within. As it stands now only rain is within the block of that if. I will fix that. It is unlikely that this causes issues but you may want to have it corrected in the borland code as well.
from crhmcode.
Thanks. I'll see if we get different values between the versions and will let Logan know about the changes.
from crhmcode.
I added some logic to consider all precipitation as rain when ratio is close to 1 and all as snow when it is close to zero.
This did have the intended effect of zeroing out snow values where previously they were very small values on the order of 10^-10 ect...
Unfortunately this did not bring evaporation, or radiation into alignment. Albedo, melftlag, and winter are not being set correctly on some steps still. I don't know yet if this is because our new calculation is not catching all near zero snow cases or if it is because my new solution alters things too much from the borland code. I am going to see if I can figure that out.
Here is the data if you want to look yourself.
Here is the altered code for reference.
double epsilon = 0.0001;
hru_snow[hh] = 0.0;
hru_rain[hh] = 0.0;
if (hru_p[hh] > 0.0) //rain or snow determined by ice bulb ratio
{
if ( abs(1.0 - ratio) <= epsilon) {
//ratio is close enough to 1 that it should 100% rain
hru_rain[hh] = hru_p[hh];
hru_snow[hh] = 0.0;
}
else if (ratio <= epsilon) {
//ratio is close enough to 0 that it should be 100% snow
hru_rain[hh] = 0.0;
hru_snow[hh] = hru_p[hh];
}
else {
//calculate rain and snow normaly
hru_rain[hh] = hru_p[hh] * ratio;
hru_snow[hh] = hru_p[hh] * (1.0 - ratio);
}
}
from crhmcode.
Thanks. I suspect that we'll still have to use a zero SWE threshold as well.
I you can push the code to your branch, I'll check it out.
from crhmcode.
I have pushed my code.
from crhmcode.
Thanks!
from crhmcode.
The albedo seems to be incorrect due to differences in the depth of SWE. This is occurring for depths >> 0.
It could be an error in the albedo decay calculations, or it could be an error in the melt rate.
It also seems that the program is not setting the Albedo to the soil value immediately when the snow depth is 0.
I'll start taking a look at the albedo and melt calculations.
OTOH, there has been some improvement in some of the previously incorrect values for net radiation as you can see for
May 1995.
When I looked at the annual discharges, some years are improved but surprisingly some are worse.
from crhmcode.
I just looked through Classabedo. It's filled with floating point tests.
This is going to take a while.
from crhmcode.
Yes, usually in cases like this it is standard practice to implement some sort of tolerance when making comparisons with floating point values.
We could do so relatively easily however the difference between single and double precision would still lead to slight differences of outcome even if we did back fill these comparisons into the borland code. Although I would guess they would be less.
However the question becomes what is the relevant tolerance for each comparison. In some cases it might make sense to have it parameterized when the tolerance represents a physical process that could be varied but if we are just trying to account for floating point issues we can just assume a small tolerance.
from crhmcode.
I have created a separate branch to explore implementing tolerances for all of the floating point comparisons. I will work on that there as I am not sure it is the best way forward but would like to explore it.
from crhmcode.
Sounds like a good idea.
Looking at the Classalbedo code, many of the comparisons are more empirically-based than physically-based.
from crhmcode.
I am noticing that this bit of code in Classalbedo.cpp line 117 is possibly causing problems:
if(net_snowD[hh] > 0.25) { // SF = SWE*2 if density 0.005 (0.5cm)
Albedo[hh] = Albedo[hh] + net_snowD[hh]*0.1*2.0; // daily value
if(Albedo[hh] > Albedo_snow[hh])
Albedo[hh] = Albedo_snow[hh];
continue;
}
There are many times where the value of net_snowD[hh] appears to be very close to 0.25 ie (0.2499999999999999) could you explain a bit about what is being tested here?
from crhmcode.
This appears to be one of the empirical relationships/tests that I was referring to.
As far as I can tell (the comment makes zero sense), it's checking the daily snowfall,
and adding a very small value to the albedo.
from crhmcode.
In the case where the daily snow is 0.2499999999999999 would it make sense to enter this branch or skip it?
from crhmcode.
I think we should be entering it
from crhmcode.
Alright, That is currently what my "fix" code is doing. However we will need to move this new style of comparison back to the borland code if we want to continue to compare the two versions.
from crhmcode.
Agreed. However, we can test the results against the old version to see if it improves things.
from crhmcode.
I did do such a test it seems that by adding a tolerance to all the comparisons in the albedo module that we moved towards closer to the borland in terms of results but I have indications that there are steps where we are moving in the opposite direction. The net effect seems to be getting closer to the borland results but there are still some places I am worried that we are "over correcting"
albedo with comparison tolerance.zip
from crhmcode.
Oof. Ok, thanks.
This module looks to be a real pain.
from crhmcode.
I don't think there is a way to get a perfect match to the borland code without accounting for these floating point differences in the borland code as well. I am going to try a few different tolerances though to see if I can get closer without that.
from crhmcode.
Thanks. It looks like there are a lot of these tests in this module
from crhmcode.
Since you are creating a new branch, it would be great to use epsilon for all of the floating-point tests in Classalbedo.
This code is critical because of the positive feedback it sets up. If the albedo is too low, then snowmelt is increased, which then further decreases the albedo.
Could we try doing this?
Thanks
from crhmcode.
That is what I did above, the data I posted is what happens when we place a comparison with tolerance in every floating point comparison in the albedo module.
from crhmcode.
Excellent. Can I try your code?
from crhmcode.
Yes just give me a moment to push to that branch.
from crhmcode.
Thanks
from crhmcode.
My code is committed on the comparison-tolerance branch. There is a logging line on common.cpp 343, and 375 that you may want to comment out before testing. Or send your output to file.
from crhmcode.
Thanks very much
from crhmcode.
I just had a meeting with Kevin Schneider and we discussed our approach to these issues dealing with floating point precision. We think it would be productive to have a conversation with you about them.
What would be a good time for you? Kevin Schneider suggested that anytime after 2:00 today or after 1:00 tomorrow would work for him.
from crhmcode.
That would be great. I have a conference call from 1-2 today, so after 2 would be good for me.
from crhmcode.
How about 2:15 then?
from crhmcode.
Works for me, thanks
from crhmcode.
Looking at Classebsm.cpp, line 205 is weird in that it is terminated with line break
double eamean = Common::estar(tmean[hh]) * rhmean[hh] / 100.0; \
This was also in the original Borland code. Any idea why it's there?
from crhmcode.
That is weird but the compiler should just ignore it as the statement is properly ended with a ';'
from crhmcode.
I did a bit more looking at the switches in Classalbedo. As we discussed yesterday at a certain point what path execution takes is a bit arbitrary but I focused in on the comparison on line 131.
Because it is a zero threshold I thought it would be particularity sensitive and correction would be easier. Also it is looking at SWE which we know often contains very small values.
After a bit of digging I found that in a 2 year period that by implementing a threshold of 1e-5 that there are 8967 instances where the SWE is above zero but under 1e-5. I wanted to see if they were evenly distributed on the HRU's and came up with this chart.
It seems HRU 5, and 10 are outliers. I don't know if that is relevant.
As you suggested above I think this really should use the same or similar threshold parameter as the pbsm module. It seems to me a place where we should change borland crhm as well.
from crhmcode.
I've been talking to John, and I think we should leave all of the thesholding alone for now.
There are too many of these things. Also, there are some issues with how we have been forcing the model
with daily precip, but hourly air temps and humiidities.
I'll send an update - an in the middle of soemthing right now.
from crhmcode.
I'm going to try re-running the model, using the old version of the code - the master branch.
This time I'm not going to use the Harder method for partitioning the precipitation into rain and snow.
I'm going to use air temperatures, which should be less susceptible to problems.
As I said above, I think that the problem is that the daily precipitation data is causing precipitation during
intervals when the RH is << 100%, which is not possible.
This is very likely the cause of snowfalls when they cannot occur.
I have run the Borland code with this configuration, so I can compare the crhmcode output against that run.
from crhmcode.
Sounds good, let me know how it goes.
from crhmcode.
I changed the precipitation partitioning to use the air temperature rather than the Harder method.
As you can see, both models now the same evaporation values:
There are some issues with other values, but I think that we can consider this issue to be closed.
from crhmcode.
Related Issues (20)
- CRHM GUI will not launch on windows HOT 24
- GUI error while running Costa et tal., 2021 prj file HOT 4
- ClassGlacier_melt_debris_cover_estimate_Mod, ClassGlacier_debris_cover, Classglacier_debris HOT 2
- Problem with installing the CLI crhm on Windows HOT 10
- segfault running CRHM with some data sets HOT 3
- glacier#5 crashes in Construct GUI
- Macro handling the parameter with Int type
- CRHMcode GUI crashes when running model with ClassWQ_pbsmSnobal module HOT 7
- Heap Corruption Detected While Running In Debug Mode in Visual Studio HOT 1
- glacier module debris-cover melt issue HOT 30
- Module hierarchy for glacier and Ayers HOT 13
- Infiltration when ZDF == 0. HOT 4
- "Manage Parameter Localy" button for module with variation HOT 1
- Variable declared as BASIN dimension not being selected correctly in CRHM_GUI Selected variable panel
- CRHM output file time stamp HOT 7
- CRHM GroupEnding name in AZ, BZ, CZ, DZ etc HOT 18
- Macro Dialog window problem for Srlab CRHM version 1.3.4 and tested version for Issue#442 HOT 2
- CRHMcode GUI - Simulation not run correctly when start date and date of first observation don't match.
- CRHM_GUI fix version for Issue_432: running crash and bug HOT 7
- Unable to compile on windows 11 in msys2 HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from crhmcode.