Incorrect psychrometric calculation causes error

[shorowit]

Issue overview

One of our models generates the following error:

   ** Severe  ** SetOutBulbTempAt: Zone Outdoor Temperatures < -100 C
   **   ~~~   ** ...check Zone Heights - Maximum Zone Height=[5].
   **  Fatal  ** Program terminates due to preceding condition(s).

The zone height and geometry are fine, so the note about checking zone heights is a red herring.

The file generates an error when I use the Dover EPW but A) the IDF works fine if I use the Wilmington EPW and B) the Dover EPW works fine with another IDF. So it’s somehow the combination of this IDF w/ the Dover EPW that generates the error.

Model: data_point.zip

Details

Some additional details for this issue (if relevant):

• Windows 10
• EnergyPlus v9.4.0

Checklist

Add to this list or remove from it as applicable. This is a simple templated set of guidelines.

• [x] Defect file added: EnergyPlusDevSupport
• [x] Ticket added to Pivotal for defect (development team task)
• [ ] Pull request created (the pull request will have additional tasks related to reviewing changes that fix this defect)

[mjwitte]

Looking at this in the debugger, the ATTIC - VENTED wetbulb temperature is the one tripping the error at -496C. This is coming from state.dataEnvrn->OutWetBulbTemp which is what the psych function in WeatherManager is returning for

OutDryBulbTemp = 1.43333
OutRelHumValue = 0.76333333333333342
OutBaroPress = 101400.00000000001

This returns OutHumRat = 0.0031902374172088472

Then Psychrometrics::PsyTwbFnTdbWPb(state, state.dataEnvrn->OutDryBulbTemp, state.dataEnvrn->OutHumRat, state.dataEnvrn->OutBaroPress); returns OutWetBulbTemp = -496 ?!

https://github.com/NREL/EnergyPlus/blob/93421b9a9d4eb1656e7fa123ded715504d107051/src/EnergyPlus/WeatherManager.cc#L1916-L1917

I think it's some anomaly in the psych function which is an iterative solution thing. Very strange. What's interesting is when I put the drybulb and dewpoint into an online psych calculator it says the wetbulb should be zero. http://daytonashrae.org/psychrometrics/psychrometrics_si.html#start

It might run if you tweak the weather data for Feb 8 hours 8 and 9 very slightly?

[amirroth]

This is one of the cached ones. Does this still happen if you disable caching? Compile with -DEP_nocache_Psychrometrics.

[mjwitte]

Yes, it happens with psychrometric caching disabled or not.

[mjwitte]

Looks like this trap may not be adequate. This is true at the time of failure, so it breaks out of the iteration loop. https://github.com/NREL/EnergyPlus/blob/93421b9a9d4eb1656e7fa123ded715504d107051/src/EnergyPlus/Psychrometrics.cc#L639-L640

[shorowit]

@JasonGlazer Any chance https://github.com/NREL/EnergyPlus/pull/8380 or https://github.com/NREL/EnergyPlus/pull/8450 would address this?

[JasonGlazer]

@shorowit I really don't know.

[jmarrec]

I traced the error by adding debug code in PsyTwbFnTdbWPb_raw in https://github.com/NREL/EnergyPlus/commit/4d350ec8c0764997b9f16b731cb48fd54adf7cb5. See output below on the timestamp where it fails:

Initial Conditions:
 TDB  =1.43333
   W  =0.00319024
Patm  =101400
tBoil =99.9948
 WBT  =1.43333

Starting iteration 1, current WBT = 1.43333
PSatstar=677.873, Patm=101400
Wstar=0.00418601
WBT >= 0,  newW = star=0.00418601, deNum = 2497.67
error = -0.000995771
WBT = 1.57667

Starting iteration 2, current WBT = 1.57667
PSatstar=684.881, Patm=101400
Wstar=0.00422958
WBT >= 0,  newW = star=0.00422958, deNum = 2497.07
error = -0.00109754
WBT = 0.0308192

Starting iteration 3, current WBT = 0.0308192
PSatstar=612.583, Patm=101400
Wstar=0.00378037
WBT >= 0,  newW = star=0.00378037, deNum = 2503.54
error = -2.2625e-05
WBT = -0.00171813

Starting iteration 4, current WBT = -0.00171813
PSatstar=611.067, Patm=101400
Wstar=0.00377096
WBT < 0,  newW = star=0.00377096, deNum = 2832.67
error = -6.75286e-05
WBT = 0.0472134

Starting iteration 5, current WBT = 0.0472134
PSatstar=613.312, Patm=101400
Wstar=0.00378491
WBT >= 0,  newW = star=0.00378491, deNum = 2503.47
error = -3.37696e-05
WBT = 0.0961602

Starting iteration 6, current WBT = 0.0961602
PSatstar=615.496, Patm=101400
Wstar=0.00379846
WBT >= 0,  newW = star=0.00379846, deNum = 2503.26
error = -6.70749e-05
WBT = -0.00241581

Starting iteration 7, current WBT = -0.00241581
PSatstar=611.032, Patm=101400
Wstar=0.00377075
WBT < 0,  newW = star=0.00377075, deNum = 2832.67
error = -6.70616e-05
WBT = -496.824
Error Trap for the Discontinuous nature of PsyPsatFnTemp function (Sat Press Curve) at ~0 Deg C

[jmarrec]

On the specific iteration 7, https://github.com/NREL/EnergyPlus/blob/5358d3eb87d1376b3484e1a38ade8274b01f0027/src/EnergyPlus/General.cc#L951-L983

EnergyPlus::General::Iterate(ResultX=0x00007fffffffaab0, Tol=0.0001, X0=-0.0024158059673152245, Y0=-0.000067061572939347645, X1=0x00007fffffffaac0, Y1=0x00007fffffffaac8, Iter=7, Cnvg=0x00007fffffffaaa0)

(lldb) p X0
(const Real64) $16 = -0.0024158059673152245

(lldb) p X1
(Real64) $17 = 0.096160198677897671

(lldb) p Y0
(const Real64) $18 = -0.000067061572939347645
(lldb) p Y1
(Real64) $19 = -0.000067074878845256852

(lldb) p std::abs(X0 - X1)
(double) $20 = 0.098576004645212892

(lldb) p DY
(Real64) $21 = 1.3305905909206939E-8

(lldb) p small
(const Real64) $23 = 1.0000000000000001E-9

(lldb) p (Y0 * X1 - Y1 * X0)
(double) $22 = -0.0000066106940700713134

[jmarrec]

Naively there's something with this line that bothers me: https://github.com/NREL/EnergyPlus/blob/5358d3eb87d1376b3484e1a38ade8274b01f0027/src/EnergyPlus/General.cc#L953

Y0 is the error. Checking for perfect equality of Y0 == 0.0 seems weird to me. The errors Y0 (current) and Y1 (previous error) are as follows when it throws:

  Y0 (error)    = -6.70616e-05
  Y1 (error-1) = -6.70749e-05

[rraustad]

I wonder if it would help to get rid of the discontinuity in curves used to calculate saturation pressure over ice and water (my guess is yes). At 0C there is a step function in PsyPsatFnTemp of about 0.06 pascals. What if instead of:

if (Tkel < 173.15) {   // Tkel < -100 C
} else if (Tkel < DataGlobalConstants::KelvinConv) {  // Tkel < 0 C
} else if (Tkel <= 473.15) { // Tkel < 200 C
)

It was:

if (Tkel < 173.15) {
} else if (Tkel < DataGlobalConstants::KelvinConvMinusOne) {
   // and in here the discontinuity is smoothed out where each endpoint of this line or curve meets the endpoints of the others
} else if (Tkel < DataGlobalConstants::KelvinConvPlusOne) {
} else if (Tkel <= 473.15) {
)

[rraustad]

At plus or minus 1 tenth C the equation is:

if (Tkel < 173.15) { } else if (Tkel < DataGlobalConstants::KelvinConvMinusOneTenth) { } else if (Tkel < DataGlobalConstants::KelvinConvPlusOneTenth) { Pascal = -12401.4579936261 + 47.6381532112526*Tkel } else if (Tkel <= 473.15) { )

[jmarrec]

That's a great idea @rraustad thanks!

[jmarrec]

I finally figured it out I think...

the ASHRAE Handbook is being slightly inaccurate in the wording, and referring to Eq 5 applying to -100°C to 0°C and Eq 6 applying to 0°C to 200°C. Quoting:

The saturation pressure over ice for the temperature range of −100 to 0°C is given by (5) The saturation pressure over liquid water for the temperature range of 0 to 200°C is given by

In fact, it is not 0°C, but the triple-point of water, which is 0.01°C. Evaluating the Eq 5 and 6 up to and from the triple-point is removing the discontinuity altogether.

[rraustad]

That's even better. Changing the transition so that there is no discontinuity should solve this.

[jmarrec]

FYI just like on #8502, the PsyPsatFnTemp function throws a warning during the iteration (the Iterate() method creates that... a bisect wouldn't). See detail of the iterations on https://github.com/NREL/EnergyPlus/issues/8502#issuecomment-864106296. I think it shouldn't do that.