HydPy-Musk-MCT calculates nan values for minimal reference discharge values #162
Description
We encountered a relevant case where HydPy-Musk-MCT can produce nan values. The following example shows the problem based on its method prepare_states:
>>> from hydpy.models.musk_mct import *
>>> simulationstep("1d")
>>> parameterstep()
>>> catchmentarea(6.25)
>>> nmbsegments(1)
>>> length(3.0)
>>> bottomslope(0.005)
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
... nmbtrapezes(1)
... bottomlevels(0.0)
... bottomwidths(0.0)
... sideslopes(2.0)
... stricklercoefficients(40.0)
>>> for e in range(30, -31, -1):
... discharge = 10**e
... sequences = model.sequences
... sequences.factors.referencewaterdepth(nan)
... model.prepare_states(discharge=discharge)
... print(e, sequences.states.courantnumber, sequences.states.reynoldsnumber)...
25 courantnumber(nan) reynoldsnumber(nan)
24 courantnumber(1.44e+19) reynoldsnumber(1.7132916434841028e+32)
23 courantnumber(1.44e+18) reynoldsnumber(1.7132916434841028e+30)
22 courantnumber(1.44e+17) reynoldsnumber(1.7132916434841028e+28)
21 courantnumber(1.44e+16) reynoldsnumber(1.713291643484103e+26)
20 courantnumber(1439999999999999.8) reynoldsnumber(1.713291643484103e+24)
19 courantnumber(144000000000000.0) reynoldsnumber(1.7132916434841027e+22)
18 courantnumber(14400000000000.0) reynoldsnumber(1.713291643484103e+20)
17 courantnumber(1440000000000.0) reynoldsnumber(1.713291643484103e+18)
16 courantnumber(144000000000.0) reynoldsnumber(1.7132916434841034e+16)
15 courantnumber(14400000000.0) reynoldsnumber(171329164348410.3)
14 courantnumber(1440000000.0) reynoldsnumber(1713291643484.103)
13 courantnumber(144000000.0) reynoldsnumber(17132916434.841028)
12 courantnumber(14400000.0) reynoldsnumber(171329164.34841028)
11 courantnumber(1440000.0) reynoldsnumber(1713291.6434841033)
10 courantnumber(144000.0) reynoldsnumber(17132.916434841027)
9 courantnumber(14400.0) reynoldsnumber(171.32916434841033)
8 courantnumber(3532.248365859715) reynoldsnumber(11.971729671764011)
7 courantnumber(1986.329226958997) reynoldsnumber(5.048436542575093)
6 courantnumber(1116.9950097524443) reynoldsnumber(2.1289080377849996)
C:\repos\HydPy\hydpy\models\wq\wq_model.py:1213: RuntimeWarning: invalid value encountered in scalar divide
fac.celerity = fac.dischargederivative / fac.surfacewidth
5 courantnumber(628.1324540151961) reynoldsnumber(0.8977530756549369)
4 courantnumber(353.2248365859712) reynoldsnumber(0.37857933294594187)
3 courantnumber(198.6329219728593) reynoldsnumber(0.15964557971452875)
2 courantnumber(111.6995005686489) reynoldsnumber(0.06732198275344577)
1 courantnumber(62.81324517287417) reynoldsnumber(0.02838944473100739)
0 courantnumber(35.32248353002075) reynoldsnumber(0.01197172957734509)
-1 courantnumber(19.863444823905166) reynoldsnumber(0.005048520551265206)
-2 courantnumber(11.170035479457479) reynoldsnumber(0.002128943296435838)
-3 courantnumber(6.289549126237608) reynoldsnumber(0.0009003019206956727)
-4 courantnumber(3.5909705454928504) reynoldsnumber(0.0003923481025973737)
-5 courantnumber(3.5754236287683545) reynoldsnumber(0.0005704997957423952)
-6 courantnumber(nan) reynoldsnumber(nan)
...
We get nan values for extremely huge and extremely tiny discharge values. The first cannot occur in reality, so we probably do not need to do anything about it. However, the latter can appear at least directly at the beginning of a model warm-up period.
The underlying problem is with the Pegasus iteration for determining the reference water depth (for tiny reference discharge values) without an initial estimate. If the reference discharge is smaller than the related tolerance, the iteration stops immediately and might return a zero reference water level for a non-zero reference discharge.
I suggest to reduce the discharge-related tolerance in case of small reference discharge values as follows:
tol_q: float = min(sol.tolerancedischarge, flu.referencedischarge[i] / 10.0)
fac.referencewaterdepth[i] = model.pegasusreferencewaterdepth.find_x(
mn, mx, 0.0, 1000.0, sol.tolerancewaterdepth, tol_q, 100
)This makes the problematic nan values disappear:
30 courantnumber(nan) reynoldsnumber(nan)
29 courantnumber(nan) reynoldsnumber(nan)
28 courantnumber(nan) reynoldsnumber(nan)
27 courantnumber(nan) reynoldsnumber(nan)
26 courantnumber(nan) reynoldsnumber(nan)
25 courantnumber(nan) reynoldsnumber(nan)
24 courantnumber(1.44e+19) reynoldsnumber(1.7132916434841028e+32)
23 courantnumber(1.44e+18) reynoldsnumber(1.7132916434841028e+30)
22 courantnumber(1.44e+17) reynoldsnumber(1.7132916434841028e+28)
21 courantnumber(1.44e+16) reynoldsnumber(1.713291643484103e+26)
20 courantnumber(1439999999999999.8) reynoldsnumber(1.713291643484103e+24)
19 courantnumber(144000000000000.0) reynoldsnumber(1.7132916434841027e+22)
18 courantnumber(14400000000000.0) reynoldsnumber(1.713291643484103e+20)
17 courantnumber(1440000000000.0) reynoldsnumber(1.713291643484103e+18)
16 courantnumber(144000000000.0) reynoldsnumber(1.7132916434841034e+16)
15 courantnumber(14400000000.0) reynoldsnumber(171329164348410.3)
14 courantnumber(1440000000.0) reynoldsnumber(1713291643484.103)
13 courantnumber(144000000.0) reynoldsnumber(17132916434.841028)
12 courantnumber(14400000.0) reynoldsnumber(171329164.34841028)
11 courantnumber(1440000.0) reynoldsnumber(1713291.6434841033)
10 courantnumber(144000.0) reynoldsnumber(17132.916434841027)
9 courantnumber(14400.0) reynoldsnumber(171.32916434841033)
8 courantnumber(3532.248365859715) reynoldsnumber(11.971729671764011)
7 courantnumber(1986.329226958997) reynoldsnumber(5.048436542575093)
6 courantnumber(1116.9950097524443) reynoldsnumber(2.1289080377849996)
5 courantnumber(628.1324540151961) reynoldsnumber(0.8977530756549369)
4 courantnumber(353.2248365859712) reynoldsnumber(0.37857933294594187)
3 courantnumber(198.6329219728593) reynoldsnumber(0.15964557971452875)
2 courantnumber(111.6995005686489) reynoldsnumber(0.06732198275344577)
1 courantnumber(62.81324517287417) reynoldsnumber(0.02838944473100739)
0 courantnumber(35.32248353002075) reynoldsnumber(0.01197172957734509)
-1 courantnumber(19.863444823905166) reynoldsnumber(0.005048520551265206)
-2 courantnumber(11.170035479457479) reynoldsnumber(0.002128943296435838)
-3 courantnumber(6.289549126237608) reynoldsnumber(0.0009003019206956727)
-4 courantnumber(3.5909705454928504) reynoldsnumber(0.0003923481025973737)
-5 courantnumber(2.0190630672110657) reynoldsnumber(0.00016540068887241715)
-6 courantnumber(1.1351995892635722) reynoldsnumber(6.972187404537214e-05)
-7 courantnumber(0.638256666401568) reynoldsnumber(2.93901978384881e-05)
-8 courantnumber(0.3554109700087458) reynoldsnumber(1.2132846193367382e-05)
-9 courantnumber(0.20322073285449246) reynoldsnumber(5.304485407223199e-06)
-10 courantnumber(0.11888762477025756) reynoldsnumber(2.4369239326058004e-06)
-11 courantnumber(0.059462051709509994) reynoldsnumber(7.971970066113762e-07)
-12 courantnumber(0.03554004537581413) reynoldsnumber(3.836496879841717e-07)
-13 courantnumber(0.020087452184504043) reynoldsnumber(1.635748077007379e-07)
-14 courantnumber(0.011366906322640714) reynoldsnumber(6.992044318061892e-08)
-15 courantnumber(0.006435270334214587) reynoldsnumber(2.991855545260473e-08)
-16 courantnumber(0.003641748267204839) reynoldsnumber(1.2790398097455229e-08)
-17 courantnumber(0.0020579481336700047) reynoldsnumber(5.451137682771535e-09)
-18 courantnumber(0.001160069122239511) reynoldsnumber(2.3107924654479305e-09)
-19 courantnumber(0.0006516366422841097) reynoldsnumber(9.721302880345536e-10)
-20 courantnumber(0.00036439077123576473) reynoldsnumber(4.049877284279415e-10)
-21 courantnumber(0.00018869721864594614) reynoldsnumber(1.428470067540782e-10)
-22 courantnumber(0.00011098223532033145) reynoldsnumber(6.638883953487566e-11)
-23 courantnumber(6.572739578666382e-05) reynoldsnumber(3.13205816352951e-11)
-24 courantnumber(3.3621302262373526e-05) reynoldsnumber(1.0757484138547101e-11)
-25 courantnumber(1.997286706636308e-05) reynoldsnumber(5.108971209242976e-12)
-26 courantnumber(1.128746069062743e-05) reynoldsnumber(2.177732039402716e-12)
-27 courantnumber(6.386762475293777e-06) reynoldsnumber(9.307239462982253e-13)
-28 courantnumber(3.6158192518855846e-06) reynoldsnumber(3.982542704697684e-13)
-29 courantnumber(2.046433689869845e-06) reynoldsnumber(1.7029720897127646e-13)
-30 courantnumber(1.1567063589899065e-06) reynoldsnumber(7.261549974031364e-14)
While working on this and searching for similar problems, I discovered that combining trapezoidal profiles and zero reference discharge can also result in nan values. Based on the above configuration:
>>> model.prepare_states(discharge=0.0)
>>> states.courantnumber
courantnumber(nan)
>>> states.reynoldsnumber
reynoldsnumber(nan)It seems fixing this can be quickly done by checking for referencedischarge == 0 instead of correctingfactor == 1.0 in Calc_CourantNumber_V1 and Calc_ReynoldsNumber_V1:
>>> model.prepare_states(discharge=0.0)
>>> states.courantnumber
courantnumber(0.0)
>>> states.reynoldsnumber
reynoldsnumber(0.0)Unfortunately, the first change has the potential to change simulation results slightly, as is the case for our no baseflow integration test. As we increase accuracy, it should tend to improve results. Additionally, it appears unlikely that the change is notable in real applications. However, it is still a change so that I will clarify this bugfix release notes (6.1.2).