Fix #8599 - Incorrect psychrometric calculation causes error #8833

jmarrec · 2021-06-21T09:27:59Z

Pull request overview

Fix Incorrect psychrometric calculation causes error #8599 - Incorrect psychrometric calculation causes error
Change PsyPsatFnTemp to correctly use the triple-point of water (0.01°C) to evaluate the two equations 5 and 6 from ASHRAE HoF thus eliminating the discontinuity

Pull Request Author

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

Reviewer

This will not be exhaustively relevant to every PR.

jmarrec · 2021-06-21T09:34:03Z

src/EnergyPlus/DataGlobalConstants.hh

-    Real64 constexpr convertJtoGJ = 1.0E-9;           // Conversion factor for J to GJ
-    int constexpr MaxSpeedLevels = 10;                // Maximum number of speed that supports
-    int constexpr ScheduleAlwaysOn = -1;              // Value when passed to schedule routines gives back 1.0 (on)
+    Real64 constexpr TriplePointOfWaterTempKelvin = 273.16; // The triple point of water, in Kelvin


New constant

jmarrec · 2021-06-21T09:34:28Z

src/EnergyPlus/General.cc

+        if (std::abs(DY) < small) {
+            DY = small;
+        }


Unrelated, but it's handy to put a debugger breakpoint... (We should have clang-format do that everywhere IMHO)

Agreed. It's standard practice but not enforced.

jmarrec · 2021-06-21T09:35:05Z

src/EnergyPlus/Psychrometrics.cc

+            if (icvg == 1) {
+                break;
+            }


Removed the "Error Trap for the Discontinuous nature of PsyPsatFnTemp function (Sat Press Curve) at ~0 Deg C." since it's no longer discontinuous.

I vaguely remember this. It looks like the unit tests hits this so good there. Gotta check it though.

Yeah, I remember this now, this is a check for where this used to have a problem. Now that the problem is corrected, this check is no longer needed.

jmarrec · 2021-06-21T09:35:18Z

src/EnergyPlus/Psychrometrics.cc

+        // Note: the ASHRAE Handbook of Fundamentals  is being slightly inaccurate in its wording,
+        // and referring to Eq 5 applying to -100°C to 0°C and Eq 6 applying to 0°C to 200°C
+        // 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.


I think you were actually the one that found this critical point. Nice!

jmarrec · 2021-06-21T09:35:36Z

src/EnergyPlus/Psychrometrics.cc

-            Real64 const C5(0.20747825e-8);                  // Coefficient for TKel < KelvinConvK
-            Real64 const C6(-0.9484024e-12);                 // Coefficient for TKel < KelvinConvK
-            Real64 const C7(4.1635019);                      // Coefficient for TKel < KelvinConvK
+        } else if (Tkel < DataGlobalConstants::TriplePointOfWaterTempKelvin) { // Tkel >= 173.15, Tkel < 273.16 (0.01°C)


The actual fix

jmarrec · 2021-06-21T09:36:13Z

src/EnergyPlus/Psychrometrics.cc

@@ -780,21 +787,21 @@ namespace Psychrometrics {

        // If below -100C,set value of Pressure corresponding to Saturation Temperature of -100C.
        if (Tkel < 173.15) {
-            Pascal = 0.0017;
+            Pascal = 0.001405102123874164;


@rraustad and I noted that the value below -100°C does not match the evaluation of the equation at -100°C. So align.

jmarrec · 2021-06-21T09:37:09Z

src/EnergyPlus/Psychrometrics.cc

+            // If above 200C, set value of Pressure corresponding to Saturation Temperature of 200C.
+        } else { // Tkel >= 173.15 // Tkel >= KelvinConv // Tkel > 473.15
+            Pascal = 1555073.745636215;
+        }


@rraustad and I noted that the value below +200°C does not match the evaluation of the equation at +200°C. So align:

One value if using the ASHRAE method: 1555073.745636215

Another value if using the IF97 method : 1554671.8682698254

(the ASHRAE and IF97 equations give slightly different values at +200°C)

jmarrec · 2021-06-21T10:01:00Z

The json big diffs: I think the epregression script is unreasonable: https://github.com/NREL/EnergyPlusRegressionTool/blob/237b8a38c280fd0ce6563b6849ae51376e326d66/epregressions/runtests.py#L630-L643

It assumes an absolute difference > 0.0001 is a big diff. This for eg:

It's categorized as a big diff, but in fact it's a 1e-12 difference in %...

In [1]: f_1 = 30310597.103115804

In [2]: f_2 = 30310597.102913737

In [3]: abs(f_1 - f_2) > 0.00001
Out[6]: True

In [7]: (f_1 - f_2) / f_1
Out[7]: 6.6665528125937005e-12

jmarrec · 2021-06-21T10:18:20Z

The RefrigeratedWarehouse one seems to be unstable to begin with. Lots of warnings about Zone Air Heat Balance is out of balance. The biggest thing I noticed in the html diff are the initialization summary diffs, and especially the number of warmup days that goes from 8 to 17. I think that may be because it would go into the error trap at iter 4 before.
Simulation time stays about the same (47s to 52s).

Most of the small mathdiffs seem to be in the EIO, Warmup Convergence Information is slightly different.

+                                                                                                                                                      v
- Warmup Convergence Information,CORE_MID,SizingPeriod: CHICAGO ANN CLG .4% CONDNS WB=>MDB,0.2742978374,0.4207863935,Pass,Pass,0.2199779749,0.5291400272,Pass,Pass
+ Warmup Convergence Information,CORE_MID,SizingPeriod: CHICAGO ANN CLG .4% CONDNS WB=>MDB,0.2742978374,0.4207863935,Pass,Pass,0.2199779749,0.5291400273,Pass,Pass

nrel-bot · 2021-07-19T00:27:21Z

@jmarrec @lgentile it has been 28 days since this pull request was last updated.

``` /home/julien/Software/Others/EnergyPlus/tst/EnergyPlus/unit/Psychrometrics.unit.cc:452: Failure The difference between result and expected_result is 496.82405055164901, which exceeds 0.001, where result evaluates to -496.82405055164901, expected_result evaluates to 0, and 0.001 evaluates to 0.001. ```

… of PsyPsatFnTemp

…rrespond to the value of the curve at these boundary temperatures

…alculated using the current method (IF97 or ASHRAE)**

…of water, not 0°C!

…p since there is none left

…nside if

…PsyPsatFnTemp change

rraustad · 2021-09-05T00:47:02Z

This is important to be included. I'd vote for sooner than later. I can't recall all the details here but I think @jmarrec has hit the nail.

rraustad · 2021-09-05T00:55:25Z

I'm OK with these changes. I am trying to think up a test to verify this works, other than a unit test. I do have a high degree of confidence in what this fixes,

rraustad · 2021-09-05T01:19:01Z

@jmarrec help me out here, it's been a while. If this corrects an issue of T=0 C (or 0 - 0.1 C), then how could I verify that this now works? I am torn between accepting this and doing some tests to make sure this is correct. I see what you did here and it looks correct, I just need to verify before sending this home. Details.

rraustad · 2021-09-05T05:30:46Z

Comparing outdoor weather data between this branch and develop doesn't show an obvious issue.

Yet, zooming in does show a difference.

rraustad

@Myoldmopar @mjwitte I recall the reason for this change and agree with these changes.

src/EnergyPlus/Psychrometrics.cc

rraustad · 2021-09-05T05:46:42Z

src/EnergyPlus/Psychrometrics.cc

            Pascal = std::exp(C1 / Tkel + C2 + Tkel * (C3 + Tkel * (C4 + Tkel * (C5 + C6 * Tkel))) + C7 * std::log(Tkel));

            // If above freezing, calculate saturation pressure over liquid water.
-        } else if (Tkel <= 473.15) { // Tkel >= 173.15 // Tkel >= KelvinConv
+        } else if (Tkel <= 473.15) { // Tkel >= 173.15 // Tkel >= TriplePointOfWaterTempKelvin


The else if above is Tkel < Real64 constexpr TriplePointOfWaterTempKelvin = 273.16; and this is <= 473.15 ... oh yeah, this is the top end so comment is wrong.

These comments are a little confusing but a second looks shows they say what they said before with updates to what was changed. I'd say this is OK.

I think the comments provide little value overall here, but there were pre-existing. and they should be correct.

rraustad · 2021-09-05T05:47:31Z

src/EnergyPlus/Psychrometrics.cc

            Pascal = std::exp(C8 / Tkel + C9 + Tkel * (C10 + Tkel * (C11 + Tkel * C12)) + C13 * std::log(Tkel));
+
+            // If above 200C, set value of Pressure corresponding to Saturation Temperature of 200C.
+        } else { // Tkel >= 173.15 // Tkel >= TriplePointOfWaterTempKelvin // Tkel > 473.15


So shouldn't this be the top end? 473.15 ?

Oh, it's only a comment issue.

Same for this comment, these are each true and true as a whole.

rraustad · 2021-09-05T05:51:18Z

Let's see what the follow up responses are before merging this one.

rraustad · 2021-09-05T05:52:48Z

I'd say this is good except for lingering comments in code.

rraustad · 2021-09-06T15:37:55Z

In the zoomed in plot above there is an outlier at X= -0.43128 (develop), Y= -0.10992 (this branch). This occurred on 4/3 19:07:30 at row 11669 (and similar point at row 11670 - 19:15:00). Plotting rows 11660 - 11680 shows the WB traces for both branches. This branch smooths out the calculated wet-bulb temperature.

Note: additional reports added to diagnose problem and corrective fix.

Output:Variable,Outside Air Inlet Node 1,System Node Humidity Ratio,detailed;
Output:Variable,Outside Air Inlet Node 1,System Node Pressure,detailed;
Output:Variable,Outside Air Inlet Node 1,System Node Wetbulb Temperature,detailed;
Output:Variable,Outside Air Inlet Node 1,System Node Dewpoint Temperature,detailed;

The outdoor conditions at this same time show that OA WB should be smooth and this branch shows the correct result. If I were to do this again I would use System Node OA temp instead of weather OAT to fill in the gaps in the zone timestep environment report. But this is good enough to show this branch corrects the problem.

rraustad

Code reviewed and tested using example file. Results show an improvement in calculated wet-bulb temperature.

jmarrec added the Defect Includes code to repair a defect in EnergyPlus label Jun 21, 2021

jmarrec self-assigned this Jun 21, 2021

jmarrec commented Jun 21, 2021

View reviewed changes

jmarrec mentioned this pull request Jun 21, 2021

JSON diff has unreasonable threshold for big diff NREL/EnergyPlusRegressionTool#70

Open

jmarrec requested review from mitchute and Myoldmopar and removed request for mitchute July 23, 2021 12:09

Myoldmopar added this to the EnergyPlus 9.6 Release milestone Aug 6, 2021

jmarrec added 13 commits August 26, 2021 07:21

Try fixing the issue #8599 - by smoothing out the discontinuity at 0C…

8b432be

… of PsyPsatFnTemp

Smooth it up by using TWO linear segment, passing by 0°C.

03f018b

Saturation pressures at the extreme (< - 100°C or > 200°C) did not co…

eb4613c

…rrespond to the value of the curve at these boundary temperatures

Set the value for T>200°C to the value that correspond to T=200°C **c…

e05a4c8

…alculated using the current method (IF97 or ASHRAE)**

Found the issue! The ashrae equations apply to/from the triple-point …

bdbaff2

…of water, not 0°C!

And now you can remove the trap for the discontinuity of PsyPsatFnTem…

3608395

…p since there is none left

Unrelated minor change: give me a spot to put a debugger breakpoint i…

62cc2cb

…nside if

Update unit test to match.

b4f7e83

clang-format [chore]

7d40fa2

Adjust unit tests that had small tolerances and were affected by the …

3118219

…PsyPsatFnTemp change

Unrelated: no reason to have these be ASSERT tests, EXPECT is fine

e1f21d6

clang format in tst [chore]

c8d9d54

jmarrec force-pushed the 8599_PsychError_SmoothPsat branch from b0fbf8f to c8d9d54 Compare August 26, 2021 05:23

Clean up after rebase with conflicts.

77bd925

jmarrec force-pushed the 8599_PsychError_SmoothPsat branch from adafb3e to 77bd925 Compare August 26, 2021 05:31

rraustad approved these changes Sep 5, 2021

View reviewed changes

rraustad reviewed Sep 5, 2021

View reviewed changes

src/EnergyPlus/Psychrometrics.cc Show resolved Hide resolved

rraustad reviewed Sep 5, 2021

View reviewed changes

rraustad approved these changes Sep 6, 2021

View reviewed changes

rraustad merged commit cf82e98 into develop Sep 6, 2021

rraustad deleted the 8599_PsychError_SmoothPsat branch September 6, 2021 16:38

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix #8599 - Incorrect psychrometric calculation causes error #8833

Fix #8599 - Incorrect psychrometric calculation causes error #8833

jmarrec commented Jun 21, 2021 •

edited by rraustad

jmarrec Jun 21, 2021

jmarrec Jun 21, 2021

rraustad Sep 5, 2021

jmarrec Jun 21, 2021

rraustad Sep 5, 2021

rraustad Sep 5, 2021

jmarrec Jun 21, 2021

rraustad Sep 5, 2021

jmarrec Jun 21, 2021

jmarrec Jun 21, 2021

jmarrec Jun 21, 2021

jmarrec Jun 21, 2021

jmarrec commented Jun 21, 2021

jmarrec commented Jun 21, 2021 •

edited

nrel-bot commented Jul 19, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021 •

edited

rraustad left a comment

rraustad Sep 5, 2021

rraustad Sep 6, 2021

jmarrec Sep 6, 2021

rraustad Sep 5, 2021

rraustad Sep 5, 2021

rraustad Sep 6, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 6, 2021 •

edited

rraustad left a comment

Fix #8599 - Incorrect psychrometric calculation causes error #8833

Fix #8599 - Incorrect psychrometric calculation causes error #8833

Conversation

jmarrec commented Jun 21, 2021 • edited by rraustad

Pull request overview

Pull Request Author

Reviewer

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jmarrec commented Jun 21, 2021

jmarrec commented Jun 21, 2021 • edited

nrel-bot commented Jul 19, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021 • edited

rraustad left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rraustad commented Sep 5, 2021

rraustad commented Sep 5, 2021

rraustad commented Sep 6, 2021 • edited

rraustad left a comment

Choose a reason for hiding this comment

jmarrec commented Jun 21, 2021 •

edited by rraustad

jmarrec commented Jun 21, 2021 •

edited

rraustad commented Sep 5, 2021 •

edited

rraustad commented Sep 6, 2021 •

edited