Tighten up functions and validity checks by dmey · Pull Request #24 · psychrometrics/psychrolib

dmey · 2019-02-24T00:57:51Z

This PR tightens up checks for handling the minimum value of humidity ratio (now reset to no lower than 1e-7) as well issues with GetTDewPointFromVapPres andGetTWetBulbFromHumRatio where two new checks MIN_ITER_COUNT and MAX_ITER_COUNT are introduced when evaluating the loop condition. MIN_ITER_COUNT checks the minimum number of iterations before exiting the loops is no less than 3. And MAX_ITER_COUNT checks that the maximum number of iterations before exiting the loops is less than 1000. These numbers can be changes.

All changes have now been propagated to all languages. All tests pass. For now I have added a simple test to check that the bounded humidity ratio is behaving as expected -- e.g.
GetTWetBulbFromHumRatio(25,1e-09,95461) == GetTWetBulbFromHumRatio(25,1e-07,95461) however do let me know if you want to add more tests and will do so.
I have already tested VBA in Excel but not committed it yet as this is dependent to my previous PR (#23). Once #23 is merged, I will commit those tests too.

@didierthevenard can you please review this and let me know for any issues. Please review #23 first. Once this is in, I can prepare the new release. Thanks

didierthevenard

I reviewed only the C version and requested some changes: not enforce a minimum number of iterations and give some kind of error if the maximum number of iterations is exceeded.

src/c/psychrolib.c

didierthevenard · 2019-03-04T02:06:24Z

src/c/psychrolib.c

-#define INVALID -99999    // Invalid value
+#define MIN_ITER_COUNT 3          // Minimum number of iterations before exiting while loops.
+
+#define MAX_ITER_COUNT 1000       // Maximum number of iterations before exiting while loops.


100 or even 20 should be sufficient. 2^20 = 1e+06 so in any bisection that we are doing we should be below the convergence threshold. NR should converge in 3 to 5 iterations given the type of very smooth function we have to work with. What's more important is that if convergence is not reached in MAX_ITER_COUNT then the program should return an error.

A couple of points about this I am not sure about:

Are you proposing to have different MAX_ITER_COUNT in the two loops?

I am not sure about throwing an error if MAX_ITER_COUNT is not reached. This is manly because there may be conditions where this is triggered but you do not necessarily want your program to stop abruptly. How about printing a warning instead?

No, one MAX_ITER_COUNT for all loops. 20 should be enough, but put 100 is you want to be extra safe. And if it does not converge in 100 iterations then there is something seriously wrong and we need to have a look, therefore it's appropriate for the program to stop or return -99999. This is just defensive programming - we never expect to get there. If we did at the moment then the program would just hang. At least with a stop or -99999 it won't.

OK for MAX_ITER_COUNT but I think that the underlying issue still remains -- for an example of no convergence you can try the following conditions:

import psychrolib as psy psy.SetUnitSystem(psy.SI) print(psy.GetTWetBulbFromRelHum(7, 0.61, 100000))

I have just looked at this in more details and I think it's now fixed in 52232e2.

The issue was with the the discontinuity of the saturated pressure curve at ~0 deg C (+/-) that leads to the calculation of dew point temperature to oscillate about the zero, thus not converging. If you assume a range of of +/-0.001, the vapour pressure at 0.001 deg C is 611.26 Pa and 611.10 Pa at -0.001 deg C I simply check if vapour pressure is within that range and break if the condition is met. Note that this will only happen when the loop index reaches the MAX_ITER_COUNT value.

psychrolib/src/python/psychrolib.py

Lines 452 to 457 in 4573c17

# The discontinuity of the saturated pressure curve at ~0 deg C (+/-) leads the

# calculation of dew point temperature to oscillate about the zero, thus not converging.

# Assuming a range of of +/-0.001, the vapour pressure at 0.001 deg C is 611.26 Pa and

# 611.10 Pa at -0.001 deg C.

if ((VapPres > 611.10) and (VapPres < 611.26) and (index == MAX_ITER_COUNT)):

break

If for some reason it had not converged after that, I throw an exception (MAX_ITER_COUNT set to 100)

psychrolib/src/python/psychrolib.py

Lines 452 to 457 in 4573c17

# The discontinuity of the saturated pressure curve at ~0 deg C (+/-) leads the

# calculation of dew point temperature to oscillate about the zero, thus not converging.

# Assuming a range of of +/-0.001, the vapour pressure at 0.001 deg C is 611.26 Pa and

# 611.10 Pa at -0.001 deg C.

if ((VapPres > 611.10) and (VapPres < 611.26) and (index == MAX_ITER_COUNT)):

break

I have now implemented this in Python only. If you are happy and there are no additional changes to make, I will propagate the changes to the other languages and add the tests to check this case too so that you can review them.

…ty-checks

didierthevenard

Checked only the c version, the issues I have in the c version are likely similar in other languages.

didierthevenard · 2019-03-16T12:36:40Z

src/c/psychrolib.c

+#define INVALID -99999            // Invalid value

-#define INVALID -99999    // Invalid value
+#define MIN_ITER_COUNT 3          // Minimum number of iterations before exiting while loops.


I thought we had agreed to drop the minimum number of iterations?

didierthevenard · 2019-03-16T12:37:10Z

src/c/psychrolib.c

-#define INVALID -99999    // Invalid value
+#define MIN_ITER_COUNT 3          // Minimum number of iterations before exiting while loops.
+
+#define MAX_ITER_COUNT 1000       // Maximum number of iterations before exiting while loops.


Should be much lower, like 100

didierthevenard · 2019-03-16T12:40:06Z

src/c/psychrolib.c

-  while (fabs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE);
+  while (((fabs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE) && (index < MAX_ITER_COUNT))
+         || (index < MAX_ITER_COUNT));
  return min(Tdp, TDryBulb);


The || (index < MAX_ITER_COUNT) part is wrong. The while loop should continue only when the first part of the test is satisfied, and the program should return an error if index > MAX_ITER_COUNT. I would actually remove the test on index from the while statement and put it in the body of the do loop.

dmey · 2019-03-16T12:53:49Z

I have just merged the changes from master but I had not acted on the minimum number of iterations yet or the other comments we had opened as I wanted to merge the fix in #26 first.
Let me add the tests first in a separate PR merge those in master, and will get back here to make the final changes. Don't review this just yet.

More changes to make -- will send an update shortly

dmey · 2019-03-16T14:49:31Z

@didierthevenard can you please look at the py version only and let me know if that's OK (It may be easier to look at the src without diff due https://github.com/psychrometrics/psychrolib/blob/dmey/tighten-validity-checks/src/python/psychrolib.py). If so, I will make the changes for the other languages. Could you also please add a code snippet for the convergence tests for VBA. Here's the py one:

psychrolib/tests/test_psychrolib_si.py

Lines 79 to 89 in bc4393c

    
           # Test that the NR in GetTDewPointFromVapPres converges. 
        
           # This test was known problem in versions of PsychroLib <= 2.0.0 
        
           def test_GetTDewPointFromVapPres_convergence(psy): 
        
               TDryBulb = np.arange(-100, 200, 1) 
        
               RelHum = np.arange(0, 1, 0.1) 
        
               Pressure = np.arange(60000, 120000, 10000) 
        
               for T in TDryBulb: 
        
                   for RH in RelHum: 
        
                       for p in Pressure: 
        
                           psy.GetTWetBulbFromRelHum(T, RH, p) 
        
               print('GetTDewPointFromVapPres converged')

didierthevenard · 2019-03-16T15:00:38Z

I had a look at the py version and I think it's great. My only recommendation is to change the
raise ValueError("Convergence not reached. Stopping.")
to
raise ValueError("Convergence not reached in GetTDewPointFromVapPres. Stopping.")
or
raise ValueError("Convergence not reached in GetTWetBulbFromHumRatio. Stopping.")
depending upon where the issue happens, so that it is easier to debug.

Will work on the vba test later today.

didierthevenard · 2019-03-16T18:54:46Z

VBA test has been added to branch djt/vba_convergence_test, I believe I have added you as a reviewer but if that has failed let me know.

…ty-checks

This PR adds test to check that the NR in GetTDewPointFromVapPres converges. This test was known problem in versions of PsychroLib <= 2.0.0

dmey · 2019-03-20T17:47:10Z

Superseded by #30. Closing.

dmey added 12 commits February 14, 2019 16:35

fortran: tighten validity checks

00a9d48

fix typo and add max iter count

295f1b0

fix fortran

9bc075c

C: tighten val checks

d335075

fix c

08e8586

js: tighten up val checks

e3dfe51

fix fortran

fc06d99

fix fortran

c11a5a4

fix logic

7df076f

VBA: tighten val checks

1f2297c

add test for clamping humratio to 10e-7

72f5f63

minor: 1E-7 -> 1e-7

155670f

dmey added the enhancement label Feb 24, 2019

dmey added this to the 2.1.0 milestone Feb 24, 2019

dmey self-assigned this Feb 24, 2019

dmey requested a review from didierthevenard February 24, 2019 00:57

dmey added 2 commits February 24, 2019 09:22

start index at 0 on all langs

7cd1bed

reduce value of MIN_ITER_COUNT: 5 -> 3.

7b09d0d

dmey added a commit that referenced this pull request Feb 24, 2019

fix notes + remove bounded humratio as changes belongs to #24

2a60bec

didierthevenard requested changes Mar 4, 2019

View reviewed changes

dmey added 4 commits March 4, 2019 08:51

remove MIN_ITER_COUNT

d24a332

fix convergence in GetTDewPointFromVapPres

52232e2

typo

4573c17

Merge remote-tracking branch 'origin/master' into dmey/tighten-validi…

e1feda8

…ty-checks

didierthevenard previously requested changes Mar 16, 2019

View reviewed changes

dmey mentioned this pull request Mar 16, 2019

GetTDewPointFromVapPres: add simple test to check convergence #27

Merged

dmey added 2 commits March 16, 2019 13:55

GetTDewPointFromVapPres: add simple test to check convergence (#27)

398c83b

add conv tests for js, py, f and c

2983b23

remove discontinuity code in saturated pressure func

bc4393c

dmey and others added 7 commits March 17, 2019 10:08

py finish adding BoundedHumRatio + modify range in tests for conv in IP

5cbae56

fix js tests by inreasing the timelimit in mocha from 2 to 40 s

e9519e2

C: integrate chnages from python version

05c5fef

F, c, js: integrate changes from python version [wip]

65fba5f

fortran make loop in GetTDewPointFromVapPres same in other impl

d56465a

Merge remote-tracking branch 'origin/master' into dmey/tighten-validi…

4f63e28

…ty-checks

Add convergence test for VBA/Excel (#28)

474e724

This PR adds test to check that the NR in GetTDewPointFromVapPres converges. This test was known problem in versions of PsychroLib <= 2.0.0

dmey closed this Mar 20, 2019

dmey deleted the dmey/tighten-validity-checks branch March 20, 2019 17:47

	# The discontinuity of the saturated pressure curve at ~0 deg C (+/-) leads the
	# calculation of dew point temperature to oscillate about the zero, thus not converging.
	# Assuming a range of of +/-0.001, the vapour pressure at 0.001 deg C is 611.26 Pa and
	# 611.10 Pa at -0.001 deg C.
	if ((VapPres > 611.10) and (VapPres < 611.26) and (index == MAX_ITER_COUNT)):
	break

Conversation

dmey commented Feb 24, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

didierthevenard left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

didierthevenard left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dmey commented Mar 16, 2019

Uh oh!

dmey commented Mar 16, 2019

Uh oh!

didierthevenard commented Mar 16, 2019

Uh oh!

didierthevenard commented Mar 16, 2019

Uh oh!

dmey commented Mar 20, 2019

Uh oh!

Reviewers

Assignees

Labels

Milestone

Development

Uh oh!

2 participants

dmey commented Feb 24, 2019 •

edited

Loading