geotopmodel / geotop

GEOtop is a distributed model of the mass and energy balance of the hydrological cycle, which is applicable to simulations in continuum in small catchments. GEOtop deals with the effects of topography on the interaction between energy balance and hydrological cycle with peculiar solutions.
GNU General Public License v3.0
38 stars 29 forks source link

Maintaining "Pi/180." instead of "FromDegToRad" #119

Closed ElisaBortoli closed 5 years ago

ElisaBortoli commented 5 years ago

PR created after noticing that three 3D tests were failing (without MeteoIO):

The executable was able to run but the results were different compared to the 2.0 version. The reason was numerical differences introduced with the definition of the new variable

constexpr double FromDegToRad = Pi/180; /** conversion from degrees to radians */

I've tried also to modify the type of the variable (i.e. `long long)` but those tests were anyway failing.

[TravisCI didn't tell me this previously since I modified the .travis.yml file to check only small_example with MeteoIO => my fault]

ElisaBortoli commented 5 years ago

@asartori86 what do you think? Did I do something wrong in the definition of the new variable? It sounds strange to me to have all those numerical errors with that substitution..

asartori86 commented 5 years ago

@ElisaBortoli, I think having a variable Pi_180 or the like to avoid to divide by 180 is recommended. It might be that the definition of FromDegToRad is inaccuare. How big are the differences? I would use M_PI from <cmath> istead of the hand written Pi

constexpr double Pi {M_PI}; // from <cmath>
constexpr double Pi_180 {Pi/180.}; // I noticed that you used `180` and not `180.` 

OR, the results are different because of the non-associativity of floating point operations

x*y/z != x*(y/z)

If the results are good or not we have to ask to the scientists

At the beginning I though of the non-associativity of floating point operations but then I looked at the faling_output file and the differences sometimes were also of order of 10^-1 or 10^-2. Precisely for the test:

Borden05m

no_reflection

I've tried also to print the values:

std::cerr << "GTConst::FromDegToRad = " << GTConst::FromDegToRad << std::endl;
std::cerr << "GTConst::Pi/180.      = " << GTConst::Pi/180. << std::endl;

but the results is the same (0.0174533).

Unfortunately no changes appear using . I think that at least for now I'll keep Pi/180 to avoid failing tests.

Thank you!