-
Notifications
You must be signed in to change notification settings - Fork 396
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Adopt the Extended Heat Index calculation in zone resilience #10548
Conversation
in TEST_F(EnergyPlusFixture, extendedHI_heatindex_compare), the HI_values used here as a baseline is output from the heatindex.py from here https://romps.berkeley.edu/papers/pubs-2020-heatindex.html
Hi @Myoldmopar and @mjwitte this is ready for review. Thanks. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, so my main concern is around the complexity of the heat index calculation. Asking functions to return pairs and tuples that include strings, and then checking the string to do different operations, is all .... a lot. I'm also highly concerned about making a copy of the entire root solver object. This needs some work. Simplification work.
I have less critical other concerns: a LaTeX error, the filename should ideally be ExtendedHI.* so that it matches our existing pattern, the unit test precision gives some feels but might be fine.
regression is not valid for extreme temperature and relative humidity conditions | ||
beyond the range of data considered by Steadman. | ||
Starting from version 24.2, the heat index calculation adopts the extended heat | ||
index method developed by Lu & Romps [17]. The previous heat index gives |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is a LaTeX error here with the & character, it needs to be escaped.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks. I just escaped it.
src/EnergyPlus/CMakeLists.txt
Outdated
@@ -262,6 +262,8 @@ set(SRC | |||
EvaporativeFluidCoolers.hh | |||
ExhaustAirSystemManager.cc | |||
ExhaustAirSystemManager.hh | |||
extendedHI.cc | |||
extendedHI.hh |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't like the lower cased names... I mean in general I am fine with it, but it looks super weird against all the other PascalCased names. I'm going to change this if I push changes up.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks. I just changed the filenames to PascalCased names.
src/EnergyPlus/extendedHI.cc
Outdated
|
||
if (show_info) { | ||
if (region == "I") { | ||
std::cout << "Region I, covering (variable phi)\n"; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are these supposed to be debug messages that can be removed? I dislike seeing std::cout directly in E+ code, even if it's only meant to be used during debugging.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
They are debugging messages. I will remove them.
src/EnergyPlus/extendedHI.cc
Outdated
// Dictionary to map eqvar_name to tuple index | ||
std::map<std::string, int> dic = {{"phi", 1}, {"Rf", 2}, {"Rs", 3}, {"Rs*", 3}, {"dTcdt", 4}}; | ||
auto eqvars = find_eqvar(state, Ta, RH); | ||
std::string eqvar_name = std::get<0>(eqvars); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are you asking a function to return a string, and then doing a string comparison to decide which code path to follow? For one thing, this arrangement seems way way overkill just to calculate heat index. And second, it cannot be a string, you need to at least use an enum. But I really feel this whole thing could be way simpler.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will simplify this
src/EnergyPlus/extendedHI.hh
Outdated
} // namespace extendedHI | ||
} // namespace EnergyPlus | ||
|
||
#endif |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No newline here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks, I just added newline at the end.
27187.715714871047, | ||
41643.766112228834, | ||
62053.2640569111, | ||
90163.72448626769}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm inclined to say all the values in this unit test are waaay too precise. We cannot know these values so precisely, so to test against them is futile. I see you have tolerance set to e-8, which is pretty tight here. But there are other ones where tolerance is more like e-4, yet the values you compare against have like 12+ decimals. It's not the end of the world, it's just a lot of numbers here :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will unify the tolerance and the values here.
These are output from this python file. I was trying to ideally match their result.
heatindex.py.zip
The e-4 part is due to the difference in the root solvers between EnergyPlus and the source python file (when EnergyPlus uses the bisection one as well).
The following is a illustration of the difference in a somewhat extreme case
Here is the difference in the code
TEST_F(EnergyPlusFixture, extendedHI_heatindex) | ||
{ | ||
|
||
// fixme: not finished |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fix me
src/EnergyPlus/extendedHI.cc
Outdated
// show_info: whether to print some messages. Might be useful in debugging | ||
// The function computes the extended heat index, in Kelvinn | ||
|
||
auto const HVACSystemRootSolverBackup = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ummmm making a copy of the whole root solver object each time we calculate heat index.... that's a no.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see. I changed it to this
auto const &HVACSystemRootSolver = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver;
auto const HVACSystemRootSolverBackup = HVACSystemRootSolver;
state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver = HVACSystemRootSolverAlgorithm::Bisection;
// do stuff
state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver = HVACSystemRootSolverBackup;
Hi @Myoldmopar, I've simplified the functions some more and addressed other comments as well. Please see if they look okay. Thanks. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is still making a copy of the root solver object in each call to heat index, and that needs to be addressed before merging. If you are needing to simply find a zero, you can do this which is what we do in a ton of places:
int constexpr maxIter = 1000;
Real64 constexpr xMin = -10.0;
Real64 constexpr xMax = 10.0;
Real64 constexpr accuracy = 0.001;
Real64 someTargetCondition = abc;
auto f = [&state] (Real64 x) { return x - someTargetCondition; };
int flag = 0;
Real64 x = 0.0;
General::SolveRoot(state, accuracy, maxIter, flag, x, f, xMin, xMax);
Is it possible to use that here rather than the HVACSolver stuff?
src/EnergyPlus/ExtendedHI.cc
Outdated
// The function computes the extended heat index, in Kelvinn | ||
|
||
auto const &HVACSystemRootSolver = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver; | ||
auto const HVACSystemRootSolverBackup = HVACSystemRootSolver; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is still making a copy of the root solver each call to heatIndex. Which happens for each zone at each time step. This is not OK. I get the feeling you shouldn't use the HVACSystemRootSolver for this. Instead just use SolveRoot?
As an aside, I still think the find_eqvar with tuple and returning different things is still a bit too complex just to calculate heat index. At this point I wouldn't hold it up just for that, but if it could be simplified further, that would be great.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi Edwin, I was using General::SolveRoot
to solve the roots. The purpose of the part involving HVACSystemRootSolver
is just so that the root searching order inside the General::SolveRoot
function follows the bisection method, which is what was used in the original python implementation. Maybe I'll just remove this part and hopefully they don't differ much when using different root finding schemes. So far the tests in the unit test file still passes if this part is removed.
I will try to simplify find_eqvar
some more
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@Myoldmopar I just changed find_eqvar
to two functions each returning a single value.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removing the
auto const &HVACSystemRootSolver = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver;
auto const HVACSystemRootSolverBackup = HVACSystemRootSolver;
state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver = HVACSystemRootSolverAlgorithm::Bisection;
...
state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver = HVACSystemRootSolverBackup;
seems to cause some NaN in the unit test HeatBalanceSurfaceManager_TestResilienceMetricReport
, when calculating this: ExtendedHI::heatindex(state, 298.15, 0.5)
. Also in unit test HeatBalanceSurfaceManager_TestThermalResilienceReportRepPeriod
the removed code was trying to set the root finding scheme, HVACRootSolverALgorithm::Bisection, following the original python implementation. This is removed to prevent copying the object
@yujiex I apologize for sending you down a distraction path here. When I saw that you were making a copy of When you said it was only 4 bytes, that of course made things different. It's just a freaking enum. Sorry for distracting you with this, please revert any unnecessary changes and let's get this merged. |
It's totally fine. It might be more clear if the variable was named something like RootSearchingMode. Currently it sounds like very memory intensive. I just reverted the changes. |
Pulling develop in here as well to get fresh CI results. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not requesting any specific changes. I left some comments that you can respond to, but I think this is probably fine to go in.
constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979 | ||
constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018 | ||
constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018 | ||
Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be constexpr-able right? And I will assume for now that these are indeed used in multiple locations where it makes sense to keep them out here in the namespace.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yes. It should be. Yeah, some stuff like Pc
is used in a few functions.
{ | ||
constexpr Real64 hc = 17.4; | ||
constexpr Real64 phi_rad = 0.85; | ||
Real64 hr = epsilon * phi_rad * sigma * (std::pow(Ts, 2) + std::pow(Ta, 2)) * (Ts + Ta); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In places where you can't constexpr
, it's still good to const
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see. Good point
Real64 m = (Pc - Pa) / (Zs(Rs) + Za); | ||
Real64 Ts; | ||
int SolFla; | ||
General::SolveRoot( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, this is the way. If you want to make it slightly clearer, you can define the lambda first so that you have plenty of room to name things nicely, and then just pass the lambda in as the argument. It's fine as an argument here though.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Make sense. Probably better as their own function in front
Tf, | ||
[&](Real64 Tf) { | ||
return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) * (Tf - Ta) / ((Zs(Rs) + Za_bar) * (Tf - Ta) + r * Ra_bar(Tf, Ta) * (Ts_bar - Tf)) - | ||
(Tc - Ts_bar) / Rs; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is definitely getting pretty long to just be a lambda argument. I'd consider making this a nicely named function with clear arguments, etc., and then just have the lambda call that function. Debugging an issue with the arithmetic in this line would be quite difficult.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, indeed. This is probably better separated out.
if (Ta == 0.0) T = 0.0; | ||
|
||
state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolver = HVACSystemRootSolverBackup; | ||
return T; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In general, this is so much better than it was, great job there. I am still quite surprised to find the sheer amount of SolveRoot calls inside this module which is just trying to calculate heat index. Anyway, thanks for the cleanups.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for reviewing. Glad I've been able to clean it out a bit. It does seem a pretty involved calculation process in their paper as well.
Rf = 2, | ||
Rs = 3, | ||
DTcdt = 4, | ||
Num |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
} | ||
} | ||
HI = (HI - 32.0) * (5.0 / 9.0); | ||
return HI; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Interesting. I guess this is fine. I start getting nervous when the amount of logic inside the unit tests grows. Like now I need to test this function as well... But I get the intent here, so I'll just move on.
constexpr Real64 Fahrenheit2Celsius(Real64 F) | ||
{ | ||
return (F - 32.0) * 5.0 / 9.0; | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I feel like we definitely have this somewhere else so you don't need your own.
377.4622562734294, | ||
381.56900731119094, | ||
385.6185252717114, | ||
389.61199074954493}}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So where did these values come from? I am guessing they are just generated from running the code, not from an external source. If so, this is OK, but not great. If these values are from a separately validated source, then that is a nice verification.
Thanks, @Myoldmopar. I will address these today |
I pushed the branch with develop pulled in and will wait for CI to give fresh results. Once that's done, this can merge. @yujiex if you want to address any of my comments via code changes, they can come later. For now just respond to them as comments however much you want. |
I see. Thanks @Myoldmopar |
|
Still happy here locally with develop pulled in, merging. Thanks @yujiex |
Thanks @Myoldmopar ! |
Pull request overview
Introduction
The heat index reflects human perceived hotness considering both temperature and humidity. The current heat index calculation in EnergyPlus is based on Steadman’s method, which gives nonphysical results in extreme weather conditions. The extended heat index fixes this issue and outputs a more realistic heat index. It can also be associated with the health outcomes (heat stroke, etc.) of an idealized human, helpful in the quantification of damages from extreme weather.
The following plot shows the absolute and percent difference between the extended HI implemented in this feature branch and the HI in the develop branch
The difference are mainly in regions with higher temperature and humidity. See the following heat index plot from the National Weather Service Website
The following is the extended heat index implemented in this feature branch.
There are also some slight differences between the extended HI implemented in this feature and the heatindex.py results by Lu and Romps. The maximum difference is around e-5. This difference is originated from the root solver difference between EnergyPlus and heatindex.py
Regression diffs
The following files have regression diffs. The ESO diffs all happen in the "ZONE ONE:Zone Heat Index [C]" variable. The table diffs are in heat index hours. Note that the heat index hours counts the number of hours with heat index in a certain range. These changes are expected as this feature changes the heat index calculation method.
Using "1ZoneUncontrolled_Win_ASH55_Thermal_Comfort.idf" as an example, the diff looks like the following
RH is 100 or close to 100. When temperature is in the range of 90F to 100F, the HI computed using the Steadman and the extended heat index method can differ by 0 to 12C.
Some other diffs are in resilience report summary tables. These tables are produced when "AllSummary" or "ThermalResilienceSummary" are present in "Output:Table:SummaryReports". For example, this is the diff of test file "_5ZoneAirCooled_LeapYear_annual.idf". The diffs are due to the change in the heat index calculation in this feature branch.
NOTE: ENHANCEMENTS MUST FOLLOW A SUBMISSION PROCESS INCLUDING A FEATURE PROPOSAL AND DESIGN DOCUMENT PRIOR TO SUBMITTING CODE
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.