Skip to content
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

Add Three Alternative Models for Clear Sky Emissivity Calculation #7562

Merged
merged 25 commits into from
Feb 19, 2020

Conversation

xuanluo113
Copy link
Contributor

@xuanluo113 xuanluo113 commented Oct 16, 2019

Pull request overview

There are two major algorithms adopted in BPS programs, one based on the work of Clark & Allen (1978), the other Martin & Berdahl (1984). Currently, EnergyPlus uses the Clark-Allen model based on a logarithmic relationship to dew point temperature, which was formulated using measurements limited to data collected over one year in San Antonio, Texas. While the original authors reported a low root mean square error (RMSE) of 10 W/m2 while used to calculate longwave radiations, studies by the International Energy Agency and Dai and Fang demonstrated that, out of the empirical models examined, Clark-Allen had among the highest high errors when tested against MODTRAN predictions and observed data, respectively. Further illustrating the limitations of Clark-Allen, a more recent assessment by Zhang et al shows that application of this model to calculate downwelling radiation in all-sky conditions tends to result in larger errors for low altitude and humid climates.

We propose to add the calibrated version of the Berdahl and Martin, Brunt and Idso models, identified as three models with higher accuracy compared to Clark-Allen and other existing popular models. The three models are calibrated using the radiation and meteorological measurements from the SURFRAD (Surface Radiation Budget Network) and ASOS (Automated Surface Observing System) operated by NOAA (National Oceanic and Atmospheric Administration). Currently seven SURFRAD stations are operating
in climatologically diverse regions over the contiguous United States including Bondville (in Illinois), Boulder (in Colorado), Desert Rock (in Nevada), Fort Peck (in Montana), Goodwin Creek (in Mississippi), Penn State University (in Pennsylvania) and Sioux Falls (in South Dakota) represent the climatological diversities.

Addresses new feature request #4999.

Pull Request Author

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

  • Title of PR should be user-synopsis style (clearly understandable in a standalone changelog context)
  • Label the PR with at least one of: Defect, Refactoring, NewFeature, Performance, and/or DoNoPublish
  • Pull requests that impact EnergyPlus code must also include unit tests to cover enhancement or defect repair
  • Author should provide a "walkthrough" of relevant code changes using a GitHub code review comment process
  • If any diffs are expected, author must demonstrate they are justified using plots and descriptions
  • If changes fix a defect, the fix should be demonstrated in plots and descriptions
  • If any defect files are updated to a more recent version, upload new versions here or on DevSupport
  • If IDD requires transition, transition source, rules, ExpandObjects, and IDFs must be updated, and add IDDChange label
  • If structural output changes, add to output rules file and add OutputChange label
  • If adding/removing any LaTeX docs or figures, update that document's CMakeLists file dependencies

Reviewer

This will not be exhaustively relevant to every PR.

  • Perform a Code Review on GitHub
  • If branch is behind develop, merge develop and build locally to check for side effects of the merge
  • If defect, verify by running develop branch and reproducing defect, then running PR and reproducing fix
  • If feature, test running new feature, try creative ways to break it
  • CI status: all green or justified
  • Check that performance is not impacted (CI Linux results include performance check)
  • Run Unit Test(s) locally
  • Check any new function arguments for performance impacts
  • Verify IDF naming conventions and styles, memos and notes and defaults
  • If new idf included, locally check the err file and other outputs

@xuanluo113 xuanluo113 added NewFeature Includes code to add a new feature to EnergyPlus IDDChange Code changes impact the IDD file (cannot be merged after IO freeze) labels Oct 16, 2019
@xuanluo113 xuanluo113 added this to the EnergyPlus Future milestone Oct 16, 2019
@xuanluo113 xuanluo113 self-assigned this Oct 16, 2019
@nealkruis
Copy link
Member

@mjwitte should this include Martin-Berdahl, too?

@xuanluo113
Copy link
Contributor Author

@nealkruis @mjwitte Great suggestion. For the first round of literature review, we selected four models, namely Berdahl & Martin, Brunt, Brutsaert and Prata as the attached document introduced in detail. These four models all most cited by literature. However, literature also suggested the newly calibrated forms of the Brunt and * Idso* model, proposed in this NFP, have better accuracy over a larger spatial domain, as they are calibrated over larger measured datasets nationwide.

If investigated in ASHRAE Standard 140, we should probably include Berdahl & Martin as well.

~Sky Emittance Models.pdf

@nrel-bot
Copy link

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

@nrel-bot-2c
Copy link

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

1 similar comment
@nrel-bot-2
Copy link

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

@xuanluo113 xuanluo113 changed the title Add Two Alternative Models for Clear Sky Emissivity Calculation Add Three Alternative Models for Clear Sky Emissivity Calculation Jan 22, 2020
Copy link
Contributor

@mjwitte mjwitte left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@xuanluo113 Some initial comments.


Alternatively, the clear sky emissivity can also be calculated using different models set by the user from options using the \textbf{WeatherProperty:SkyTemperature} object, including the calibrated version of the Berdahl and Martin, Brunt and Idso models (2017), listed as follows.

Calibrated Berdahl and Martin model:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor, but I've usually seen this referenced as Martin Berdahl rather than Berdahl Martin. Please check recent literature to see what's the most common citation of this model. It may depend on whether one is referencing this 1983 LBNL report https://escholarship.org/content/qt71v36429/qt71v36429.pdf or the 1984 paper you've referenced. Are the equations the same in both references?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

They use the same equation for sky emissivity.

Emissivity of clear skies - Berdahl and Martin (https://doi.org/10.1016/0038-092X(84)90144-0) (cited 142) is the first to introduced the emissivity equation. Martin and Berdahl (https://doi.org/10.1016/0038-092X(84)90162-2) (cited 158) introduced the sky radiant temperature equation from this emissivity estimation.

I still tend to call it Berdahl and Martin because ESP-r that currently uses this algorithm calls it Berdahl and Martin algorithm. (http://www.esru.strath.ac.uk/Programs/ESP-r_FAQ.htm)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok

idd/Energy+.idd.in Show resolved Hide resolved
src/EnergyPlus/WeatherManager.cc Show resolved Hide resolved
@@ -839,7 +839,7 @@ \subsubsection{Inputs}\label{inputs-7-013}

\subsection{WeatherProperty:SkyTemperature}\label{weatherpropertyskytemperature}

Sky Temperature, or radiative sky temperature, is internally calculated by EnergyPlus using an algorithm using horizontal infrared radiation from sky, cloudiness factors and current temperature. The algorithm is fully described in the Engineering Reference document. For flexibility, the following object can be entered to override the internal calculations. Much of the literature describes the sky temperature as relative to either drybulb or dewpoint temperature.
Sky Temperature, or radiative sky temperature, is internally calculated by EnergyPlus using an algorithm using (1) horizontal infrared radiation or (2) an empirical model using sky cloudiness factors and current clear sky emissivity. By default, EnergyPlus calculates clear sky emissivity using Clark-Allen model, and its algorithm is fully described in the Engineering Reference document. For flexibility, the following object can be entered to adopt alternative sky emissivity calculation methods, or to override the entire internal sky temperature calculation from schedule import. Alternative methods of sky emissivity calculation include the calibrated forms of Martin \& Berdahl, Brunt, and Idso model, and their algorithms are described in the Engineering Reference as well. Much of the literature describes the sky temperature as relative to water vapor pressure, drybulb or dewpoint temperature.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This sentence reads very differently from the original regarding horizontal infrared radiation. The original implies HorizIR is used in conjunction with the empirical model, but the new text implies HorizIR is a separate method (which doesn't have an option) so I'm confused by the revision.

Oh, I see from the code that it will use horizontal IR if present in the weather data, if one of the models is selected, but not if it's scheduled. That needs to be much more clear here. And should there be an option to always use the model even when horizontal IR is present? Not sure which way I lean on that question.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The logic is re-illustrated.


N = opaque sky cover \{tenths\}
where, N = opaque sky cover \{tenths\}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Applicable references from the NFP need to be added to the Engineering Ref section.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

SkyTempScheduled = true;
}

if (!SkyTempScheduled) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why all the extra gymnastics above? Seems that lines 3462:3466 could be eliminated and this line could be
if (!Environment(Envrn).WP_Type1==WP_ScheduleValue)
Oh, but I see WP_Type1 is the index of the WeatherProperty. So, seems like it would be more efficient to change WP_Type1 to be WP_SkyTempType and set it when processing the WeatherProperty object? Then all of these ifs could be more direct.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a SkyTempModel field in the Environment. It's still required in WeatherProperty object for matching with Environments and reading in the schedules.

src/EnergyPlus/WeatherManager.cc Outdated Show resolved Hide resolved
@@ -107,6 +107,9 @@ namespace WeatherManager {
extern int const WP_ScheduleValue; // User entered Schedule value for Weather Property
extern int const WP_DryBulbDelta; // User entered DryBulb difference Schedule value for Weather Property
extern int const WP_DewPointDelta; // User entered Dewpoint difference Schedule value for Weather Property
extern int const WP_BruntModel; // Use Brunt model for sky emissivity calculation
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do these all need to be extern? Looks like they're only used in WeatherManager.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mjwitte I "extern" them for the unit test only (otherwise they cannot be declared in the unit test file). Any better suggestions?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The ideal option would be to just use an enum class for options like this. There are occasions where that is not feasible, such as if a variable of that enum type is used as an integer output variable. But those are very few, and in most cases, changing to an enum can fix this. I'm not pushing for that now, but in the future, that would be great.

WPSkyTemperature(Item).IsSchedule = true;
WPSkyTemperature(Item).SchedulePtr = Found;
}
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not finding any code that reports the selected sky temperature model in the eio output. That should have happenend when WeatherProperty:* was added. Since this can be set per environment, it could be added to this eio output:

! <Environment>,Environment Name,Environment Type, Start Date, End Date, Start DayOfWeek, Duration {#days}, Source:Start DayOfWeek,  Use Daylight Saving, Use Holidays, Apply Weekend Holiday Rule,  Use Rain Values, Use Snow Values

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

@mjwitte
Copy link
Contributor

mjwitte commented Feb 7, 2020

@xuanluo113 Thanks for addressing comments. Is this ready for final review?

@xuanluo113
Copy link
Contributor Author

@xuanluo113 Thanks for addressing comments. Is this ready for final review?

@mjwitte Yes. Thanks, Mike.

@Myoldmopar
Copy link
Member

EIO changes look to be in order now. @mjwitte can you easily rerun the experiment to ensure a default weather property object gives the same results? If not I am happy to do it. CI otherwise looks happy. I'm going to look over the code.

Copy link
Member

@Myoldmopar Myoldmopar left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had some comments, but none urgent enough to necessitate pushing more commits unless there is other stuff to commit.

@@ -107,6 +107,9 @@ namespace WeatherManager {
extern int const WP_ScheduleValue; // User entered Schedule value for Weather Property
extern int const WP_DryBulbDelta; // User entered DryBulb difference Schedule value for Weather Property
extern int const WP_DewPointDelta; // User entered Dewpoint difference Schedule value for Weather Property
extern int const WP_BruntModel; // Use Brunt model for sky emissivity calculation
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The ideal option would be to just use an enum class for options like this. There are occasions where that is not feasible, such as if a variable of that enum type is used as an integer output variable. But those are very few, and in most cases, changing to an enum can fix this. I'm not pushing for that now, but in the future, that would be great.

tst/EnergyPlus/unit/WeatherManager.unit.cc Outdated Show resolved Hide resolved
"Dewpoint Difference Schedule Value",
"Brunt",
"Idso",
"Berdahl and Martin"});
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would have preferred a vector here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Myoldmopar like a std::vector? I thought this should be a const array.

double TDewC = min(DryBulb, DewPoint);
double SatPress = PsyPsatFnTemp(DryBulb) * 0.01;
double PartialPress = RelHum * SatPress;
double TDewK = TDewC + TKelvin;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Marking these as const would be best-practice

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why assign all of these here when only some are used depending on the sky model type? Seems these should all be moved within the respective sky model where they are used.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's also true @mjwitte. Though PartialPress (and therefore SatPress) is used by two of them, and TDewC is implicitly used in two of them, so there would be some duplication. I'm OK with leaving it like this or putting the initialization into each IF block.

src/EnergyPlus/WeatherManager.cc Show resolved Hide resolved
@mjwitte
Copy link
Contributor

mjwitte commented Feb 13, 2020

@Myoldmopar I'll retest today/tonight.

@mjwitte
Copy link
Contributor

mjwitte commented Feb 13, 2020

@xuanluo113 @Myoldmopar 1ZoneUncontrolled_DDChanges still crashes here in debug build on Windows. Will debug here an let you know what I find.

@xuanluo113
Copy link
Contributor Author

1Zone.zip

Thanks, Mike and Edwin. Attached the local test cases and weather files for reference.

src/EnergyPlus/WeatherManager.cc Outdated Show resolved Hide resolved
double TDewC = min(DryBulb, DewPoint);
double SatPress = PsyPsatFnTemp(DryBulb) * 0.01;
double PartialPress = RelHum * SatPress;
double TDewK = TDewC + TKelvin;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why assign all of these here when only some are used depending on the sky model type? Seems these should all be moved within the respective sky model where they are used.

@@ -3731,6 +3750,27 @@ namespace WeatherManager {
return (fmod(interpAng, 360.)); // fmod is float modulus function
}

Real64 CalcSkyEmissivity(int ESkyCalcType, Real64 OSky, Real64 DryBulb, Real64 DewPoint, Real64 RelHum){
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why have RelHum as an argument when DryBulb and DewPoint are also here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RelHum is used to calculate the partial water pressure directly. Otherwise, kinda have to convert it inside function.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@xuanluo113 you'll need to address the crash and the question from @mjwitte about the function signature.

CalcSkyEmissivity is declared as:

Real64 CalcSkyEmissivity(..., Real64 DryBulb, Real64 DewPoint, Real64 RelHum);

With relative humidity as the last argument. You also have dry bulb and dew point so unless you are saying the barometric pressure is changing, you only need to pass in two psychrometric values to locate your state. However, even with that variable there, you are calling it with humidity ratio:

ESky = CalcSkyEmissivity(..., DryBulb, TomorrowOutDewPointTemp(TS, Hour), OutHumRat);

And that's causing a crash because OutHumRat is uninitialized at that point.

Please address that as soon as possible and make sure to pull develop in right before you push your changes so we get clean CI results.

@Myoldmopar @mjwitte Yes. I'm working on it now. But the calculation of partial vapor pressure requires relative humidity directly, while the other places requires dew point. If I only pass two psychrometric vars, I have to re-do the psychrometric conversion again inside function, which is a kind of duplicate calculation?

@mjwitte
Copy link
Contributor

mjwitte commented Feb 13, 2020

Confirmed that 5ZoneAutoDXVAV and 5ZoneAutoDXVAV with a default WeatherProperty:SkyTemperature object produced identical results.

@Myoldmopar
Copy link
Member

@xuanluo113 you'll need to address the crash and the question from @mjwitte about the function signature.

CalcSkyEmissivity is declared as:

Real64 CalcSkyEmissivity(..., Real64 DryBulb, Real64 DewPoint, Real64 RelHum);

With relative humidity as the last argument. You also have dry bulb and dew point so unless you are saying the barometric pressure is changing, you only need to pass in two psychrometric values to locate your state. However, even with that variable there, you are calling it with humidity ratio:

ESky = CalcSkyEmissivity(..., DryBulb, TomorrowOutDewPointTemp(TS, Hour), OutHumRat);

And that's causing a crash because OutHumRat is uninitialized at that point.

Please address that as soon as possible and make sure to pull develop in right before you push your changes so we get clean CI results.

@xuanluo113
Copy link
Contributor Author

xuanluo113 commented Feb 14, 2020

@mjwitte @Myoldmopar I fixed the urgent ones and will wait for the CI. Three remaining issues:

  1. RelHum is calculated in the caller. I'm not sure if it's better to (1) recalculate it again in the callee when used or (2) passing an extra argument sometimes not used.
  2. I thought const array of string here could perform better than vector of strings.
  3. Mike suggested adding another field in IDD Use Weather File Horizontal IR, which is implicitly yes right now. I only included the IDD change in this commit and is working on it while waiting for CI.

@mjwitte
Copy link
Contributor

mjwitte commented Feb 14, 2020

  1. It's fine to pass RelHum rather than recalculate it.
  2. I defer to @Myoldmopar
  3. Thanks.

@Myoldmopar
Copy link
Member

  1. Yeah, leave RelHum in there as well
  2. I'm saying making it a const vector of strings. We are moving away from Array1D soon, so adding new instances is counterproductive, that's all. Honestly though, at this point, just leave it as-is.
  3. 👍

@xuanluo113
Copy link
Contributor Author

@mjwitte @Myoldmopar The CI is clear and the eio is showing.

The small eso diff is due to the line I added the here, which I think should be added.

@mjwitte
Copy link
Contributor

mjwitte commented Feb 17, 2020

@xuanluo113 I think I'm ok with the change to calculate TomorrowOutRelHum in the else block, but I'm puzzled trying to find where it was set previously such that the eso output in develop is rounded to a single decimal place, but this calculation is not. For example, what was 45.0 is now 44.9999999999.

Copy link
Contributor

@mjwitte mjwitte left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tested various options (ClarkAllen, default, Brunt, ClarkAllen and Brunt with epw HorizIR=No), all looks good, saw changes in sky temp where expected, eio outputs were consistent with the input choices.

We're almost there.

  1. Please add technical reference comments in WeatherManager::CalcSkyEmissivity
  2. Add an entry for this in \src\Transition\OutputRulesFiles\OutputChanges9-2-0-to-9-3-0.md describing the new field at the end of the eio "Environment" lines.
  3. The engineering ref needs some cleanup, exponent at the end of eq 5.5, and the subscript on the left side of eq 5.6.
    image

Copy link
Contributor

@mjwitte mjwitte left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Confirmed requested changes. This is good to go.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
IDDChange Code changes impact the IDD file (cannot be merged after IO freeze) NewFeature Includes code to add a new feature to EnergyPlus
Projects
None yet
Development

Successfully merging this pull request may close these issues.

9 participants