Fix overlapping shadow calculations for shading surface transmittance >0 and <1 #10040

mjwitte · 2023-05-23T22:05:52Z

Pull request overview

Fixes Fractional shading surface transmittance schedule gives strange results (but not always) #10022
Corrects the transmittance calculation for overlapping shadow areas when one or both of the shadows has a transmittance between 0 and 1.
This changes the solar gains on surfaces that are shaded by partially transmitting surface(s) with overlapping shadows.
There is no change for times when the shading surface transmittance is either 0 (opaque) or 1 (transparent).
Some code cleanup in SolarShading::DeterminePolygonOverlap.
Added some missing initial values in SolarShading.hh.
Updated the engineering reference section on overlapping shadows.

Diffs

Diffs are expected for all test files based on SolarShadingTest and SurfaceTest, because these have shading surfaces with a transmittance schedule.

eio diffs occur even in the files that use imported surface shading, because the sky diffuse view factors reported in the eio are impacted by the fix.

The theory

When overlapping shadows impact a surface, the surface sunlit area is calculated by subtracting each shadow, and then adding back the area of each overlap to avoid subtracting the area twice. This works great for opaque shadows, but for partially transmitting shadows, the transmittance of the overlap was the product Tau1*Tau2. That sounds correct, because that's what you get with two overlapping shades, of say 50% and 20% transmittance, you should get 10% after passing through both. But that's what the final result should be, not the value for the correcting term in the overlap area. Here's a screen shot from a spreadsheet that shows derivation of the new formula and some check calculations (attached if you want to play with values: 10022 Shading Frac Details.xlsx).

The updated Engineering Reference section "Overlapping Shadows" offers another explanation of this. Changes begin after Figure 5.11: Complex Overlapping Condition.

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.

mjwitte · 2023-05-24T17:30:35Z

Looking back in the archives, this same issue was addressed long ago.
Fractional transmittance with multiple shading devices may not be accurate. (CR #7383) #2167

Found 2 relevant commits.
https://github.com/NREL/EnergyPlusArchive/commit/1df5c16892793cc33eaf14018fc34fa8593da085
This added some logic in DeterminePolygonOverlap to use (1-HCT3) in some cases:

      HCT(NS3) = HCT(NS2)*HCT(NS1)  ! Determine transmission of overlap
+        if (HCT(NS2) /= 1.0 .and. HCT(NS2) /= 0.0 .and. HCT(NS1) /= 1.0 .and. HCT(NS1) /= 0.0) THEN
+          HCT(NS3)=1.0-HCT(NS3)

https://github.com/NREL/EnergyPlusArchive/commit/56e08a151128f55417e888b47dbf6b0ed4a9af29 added some more logic:

      HCT(NS3) = HCT(NS2)*HCT(NS1)  ! Determine transmission of overlap
        if (HCT(NS2) /= 1.0 .and. HCT(NS2) /= 0.0 .and. HCT(NS1) /= 1.0 .and. HCT(NS1) /= 0.0) THEN
+          if (HCT(NS2) >= .5 .and. HCT(NS1) >= .5) then
            HCT(NS3)=1.0-HCT(NS3)

This almost matches what's currently in DeterminePolygonOverlap

            Real64 const HCT_1(state.dataSolarShading->HCT(NS1));
            Real64 const HCT_2(state.dataSolarShading->HCT(NS2));
            Real64 HCT_3(HCT_2 * HCT_1); // Determine transmission of overlap
            if (HCT_2 >= 0.5 && HCT_1 >= 0.5) {
                if (HCT_2 != 1.0 && HCT_1 != 1.0) {
                    HCT_3 = 1.0 - HCT_3;
                }

The logic was rearranged slightly in da5f9fd#diff-18fbb8a5b1334b62513de2441bdb4bb4ef1fbe5fbec24d5dd0effcd019dfc33dL3216.

So, the question becomes, what's wrong with this logic?

…ngSched

…unit tests

rraustad · 2023-06-04T03:47:22Z

src/EnergyPlus/SolarShading.cc

+            Real64 const HCT_1 = state.dataSolarShading->HCT(NS1);
+            Real64 const HCT_2 = state.dataSolarShading->HCT(NS2);
+            Real64 HCT_3 = HCT_1 * HCT_2;
+            if (HCT_2 != 1.0 && HCT_1 != 1.0) {


Isn't this the same as if (HCT_3 != 1.0)

Hmm, at HCT_1 and HCT_2 = 1, HCT_3 = (HCT_1 * HCT_2) = 1 and HCT_3 = (HCT_1 + HCT_2) - (HCT_1 * HCT_2) = 1. And at all other transmittances, HCT_3 = (HCT_1 + HCT_2) - (HCT_1 * HCT_2). Why the need for transmittance checking?

I compared this branch to this and there are diffs. I assume it is because either HCT_1 or HCT_2 is 1 but not both where (1 + x) - (1 * x) = 1 and would have been x before. If HCT_1 and HCT_2 are 1 and 0.5, the answer is 0.5. If 0.999999 and 0.5 are used, the answer is 0.9999995. Why is that?

Real64 const HCT_1 = state.dataSolarShading->HCT(NS1); Real64 const HCT_2 = state.dataSolarShading->HCT(NS2); state.dataSolarShading->HCT(NS3) = (HCT_1 + HCT_2) - (HCT_1 * HCT_2);

This same function gets used for a variety of cases including: two overlapping shadows on an exterior surface, and shadows that penetrate a window onto an interior surface (for full interior solar dist), and some special overlaps that force once of the HCT_* values to 1.0 in the function parameters which I think is intended to just calculate the overlap area and not care about the transmittance.

Real64 HCT_3 = HCT_1 * HCT_2; if (HCT_2 != 1.0 && HCT_1 != 1.0) { HCT_3 = (HCT_1 + HCT_2) - HCT_3; }

The overlap transmittance calc always needs the product of HCT_1 * HCT_2.
If either HTC_2 == 1.0 OR HCT_3 == 1.0, then we're done and skip the more complex calc.
But if they are both != 1.0 then we're doing overlapping shadows and we need the second expression. I originally tried just doing the full expression without checking for 1.0 and it introduced unwanted diffs that I can't fully explain.

I suppose zero is another special case, if either HCT_1 or HCT_2 is zero, then HCT_3 = the nonzero one
So maybe there's a more efficient way to code this. This gets called frequently, so it's worth optimizing.

…ngSched

mjwitte · 2023-06-07T16:16:37Z

Here's the before and after using the defect files 1ZoneUncontrolled_Win_AddOverhang and 1ZoneUncontrolled_Win_AddOverhangAndSouthShadeScheduled.

For the base wall on the winter day, with develop the sunlit area flattens out mid-day at 0.8 (the south shade transmittance value), ignoring the impact of the overhang. With this branch, the impact of both overhang and partially transmitting shade are accounted for.

For the window on the winter day, the overhang has a tiny impact, so develop and this branch both show the correct result of 0.8 from the south shade.

For both the wall and the window, on the summer day, with develop adding the scheduled south shade increases the sunlit fraction instead of decreasing it. This branch shows the impact of both.

mjwitte

Here's a walk-thru. I'll merge in (fast-moving) develop once again and push for final reviews.
Requesting reviews from @rraustad, @Myoldmopar, @vidanovic, and @nealkruis

mjwitte · 2023-06-07T16:50:44Z

src/EnergyPlus/SolarShading.cc

+    // Check for exceeding array limits.
    if (NS3 > state.dataSolarShading->MaxHCS) {


Note this if at the top of the function which throws a warning and returns. There were two places further down that checked the same condition (or the opposite) which have been removed now, because they were moot or even dead.

mjwitte · 2023-06-07T16:51:40Z

src/EnergyPlus/SolarShading.cc

-    if (NV3 < state.dataSolarShading->MaxHCV && NS3 <= state.dataSolarShading->MaxHCS) {
+    if (NV3 < state.dataSolarShading->MaxHCV) {


If the second condition was false, we would have returned already up above.

This also looks correct.

mjwitte · 2023-06-07T16:52:21Z

src/EnergyPlus/SolarShading.cc

@@ -4791,26 +4785,6 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
            state.dataSolarShading->TrackTooManyVertices(state.dataSolarShading->NumTooManyVertices).SurfIndex2 =
                state.dataSolarShading->CurrentSurfaceBeingShadowed;
        }
-
-    } else if (NS3 > state.dataSolarShading->MaxHCS) {


This entire else if duplicates the block at the top of the function. This is dead code.

NS3 is const so that didn't change. As long as MaxHCS hasn't changed (from a call in this function like ClipPoly?) then I agree.

MaxHCS is a user input in ShadowCalculation "Maximum Figures in Shadow Overlap Calculations". It's never changed, just used to set some array sizes and then checked against in multiple places.

mjwitte · 2023-06-07T16:53:17Z

src/EnergyPlus/SolarShading.cc

-            if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0)
+            if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0) {


Just changing a one-line if to multi-line for clarity.

mjwitte · 2023-06-07T16:57:19Z

src/EnergyPlus/SolarShading.cc

-    } else if (NV3 > state.dataSolarShading->MaxHCV) {
+    } else {


This if is checking the number of vertices (not MaxHCS) but a simple else is adequate here.

mjwitte · 2023-06-07T16:57:47Z

src/EnergyPlus/SolarShading.cc

+            if (HCT_2 == 1.0 || HCT_1 == 1.0) {
+                state.dataSolarShading->HCT(NS3) = HCT_1 * HCT_2;
+            } else {
+                // Determine transmission of overlap which corrects for prior shadows
+                // The resulting transmission of overlapping shadows is HCT_1 * HCT_2
+                // Shadows with HCT_1 and HCT_2 have already been applied in the overlapping area
+                // so the correction is the difference between (HCT_1+HCT_2) and (HCT_1*HCT_2)
+                state.dataSolarShading->HCT(NS3) = (HCT_1 + HCT_2) - HCT_1 * HCT_2;


Here's the heart of the fix.

mjwitte · 2023-06-07T16:59:15Z

src/EnergyPlus/SolarShading.hh

-    int FINSHC;                     // HC location of first back surface overlap
-    int FRVLHC;                     // HC location of first reveal surface
-    int FSBSHC;                     // HC location of first subsurface
+    int FBKSHC = 0;                 // HC location of first back surface


Added some missing initial values that were the root cause of differing unit test results on linux and windows. The unit tests were also updated.

mjwitte · 2023-06-07T17:01:50Z

tst/EnergyPlus/unit/SolarShading.unit.cc

@@ -3973,3 +3974,1055 @@ TEST_F(EnergyPlusFixture, ShadowCalculation_CSV)
        EXPECT_EQ(expected_values, stream_str);
    }
 }
+
+TEST_F(EnergyPlusFixture, SolarShadingTest_PolygonOverlap)


Three new unit tests.
SolarShadingTest_PolygonOverlap - tests similar to the 1Zone defect file with an overhang and a tree.
SolarShadingTest_PolygonOverlap2 - tests two large south shades that both fully shade the wall and window, one behind the other. The expected results are simply the product of the transmittances.
SolarShadingTest_PolygonOverlap3 - tests three large south shades, just to be sure the method works for more overlaps.

In case you're interested, here's the test failures with develop:
PR10040-ShadingTransmittance-UnitTestFailureswithDevelop.txt

…ngSched

rraustad · 2023-06-07T19:01:18Z

doc/engineering-reference/src/climate-sky-and-solar-shading-calculations/shading-module.tex

+\end{equation}
+
+The actual $\tau$ of overlapping polygons i and j is the product of $\tau$\(_{i}\) and $\tau$\(_{j}\). For example, if a wall is fully shaded by two partially transparent shades, one with a transmittance of 0.8 and one with a transmittance of 0.5, the resulting sunlit fraction would be $0.8 * 0.5 = 0.4$. Because each shadow is first applied with its own transmittance, the transmittance for the overlap in the summation is the difference between the sum and the product of the shadow transmittances:
+


This is the only issue I have with these changes. And is the fact that I can't resolve this math. The first sentance gives the correct answer of 0.4 and I agree with that transmittance calculation. The second sentance defines that calculation as (0.8+0.5)-(0.8*0.5)=0.9. The figure you show in the PR discussion looks correct but I can't get past this math. I guess somewhere after this calculation there is more calculatioins that use the original (0.8+0.5)=1.3 which is then reduced by the result here of 0.9 giving 0.4 as the final answer. Could you show that same figure using just (HCT_1*HCT_2) as the result (or add it to the figure)? At least then, if the result looks incorrect, I could sleep. I did look at the unit test and need to look again for the sunlit area fraction of the "overlap area".

@rraustad Valid questions. See the "The theory" here. If those example calculations are helpful, I can add that to the engineering ref. If not, I'll try another way to illustrate it.

You can't argue with a spreadsheet! In the derivation spreadsheet, seeing a "transmittance" of 0.88 when it is obviously 0.32 is mind blowing. An area is not being calculated, nor is a transmittance. The sunlit area fraction is being calculated. So:

What you have already said in the docs, plus: "The overlapped area (transmittance?) correction factor used to calculate sunlit area fraction is derived from the individual surface transmittance values as (your formula)".

or, "When shading surfaces overlap, the sunlit area fraction is calculated using the shading surface transmittances and a correction for the overlapping shaded area. The area weighted sunlit area fraction is calculated as the unshaded area, plus shading area 1 multiplied by shading 1 surface transmittance, plus shading area 2 multiplied by shading 2 surface transmittance (note here the overlapped area has been included twice), plus the overlapped area multiplied by a correction factor of ()." you get the idea. Whatever you like, just don't call the overlap multiplier a transmittance.

rraustad · 2023-06-08T10:54:55Z

doc/engineering-reference/src/climate-sky-and-solar-shading-calculations/shading-module.tex

-\begin{longtable}[c]{@{}ll@{}}
-\caption{Surface / Area Characteristic / Convention \label{table:surface-area-characteristic-convention}} \tabularnewline
+Each shadow's area is subtracted from the receiving surface and so on through multiple overlaps where the sign of the overlap area is the product of the signs of the overlapping areas. For the shadows in Figure~\ref{fig:multiple-shadow-overlaps}, start with the receiving surface area A, subtract shadow area B, then subtract shadow area C. Because shadows B and C overlap, the overlap area D has been subtracted twice, so add back area D to get the final sunlit area on the receiving surface.
+


I understand the math now. The key is the sign convention when adding in the area * transmittance values. The 2 shadows "area * trans" are "subtracted" from the total area. That winds up being a (1 - HCT_1 + HCT_2) term for the overlap. Then the overlap is "added" back to that value using the difference of applied trans - actual trans which gets back to actual trans as (1 - HCT_1 + HCT_2) + (HCT_1 + HCT_2) - (HCT_1 * HCT_2) for the overlap. This yields (1 - (HCT_1 * HCT_2)) as the sunlit fraction for the overlap area. This is the original equation in the model. But that is not what needs to get added back as @mjwitte has shown. The overlap area needs to add back (HCT_1 + HCT_2) - (HCT_1 * HCT_2) to get back to the correct value which @mjwitte has said multiple times. This model therefore calculates the resulting sunlit area fraction based on shading surface transmittances.

As another test, I used 2 shading surfaces with 0.7 transmittance that covered the entire area and had a miniscule overlap. I would expect a sunlit area fraction of 0.7. The spreadsheet shows 0.699998.

@rraustad Welcome to my rabbit hole! I'm working on some doc additions, including sample calculations a la spreadsheet.

I did make sure to check the other extreme here with overlapped shades with area 19.9999 with miniscule NON overlap (D = 19.9998) and expected a result of 0.49. The actual result was 0.490002.

rraustad

Everything in the branch makes sense. Whether doc updates are needed or not is left to the developer.

vidanovic · 2023-06-08T14:53:20Z

@mjwitte I do not see any changes in sun patches coming through the windows which is good.

vidanovic

Suppose for my window tests I should add a file that will have exterior overhang so that incoming rays are tested as well. As already said, I do not see changes at the interior patch distributions and that is good.

…ngSched

rraustad · 2023-06-09T18:00:26Z

These are the diffs from SolarShading eio reports for Heat Transfer Surfaces. These type of diffs are expected with the change to the overlap model sunlit area calculation.

rraustad · 2023-06-09T18:06:34Z

Unit tests run successfully. @Myoldmopar this is ready to merge.

Myoldmopar · 2023-06-14T14:44:27Z

I'm just doing a quick local build and test of this, and if it's all clean I'll merge it next.

Myoldmopar · 2023-06-14T14:59:47Z

Everything is happy, merging. Thanks @mjwitte .

As an aside, this set of CI results on develop should be back to all clean except the SSC Battery failure on the debug build, which may be commented out soon because it's not helpful just always failing.

ecoeficiente · 2023-07-31T08:54:49Z

Hi, does this also fix this?
Solar through multiple partially transmitting shading surfaces is not correct at the receiving surface #6420

mjwitte · 2023-07-31T19:28:54Z

Hi, does this also fix this? Solar through multiple partially transmitting shading surfaces is not correct at the receiving surface #6420

Yes, it does! Here's plot of the results for the roof shaded just by one flat shade at 0.1 transmittance, compared to the roof with the 0.1 flat shade plus a 0.1 transmittance box above it, with an expected shading fraction of 0.01 when the sun is overhead (the shading surfaces are slightly above the roof, so, some solar gets in from the side). You can see how v23.1 was clearly wrong (more solar with the added shades), and v23.2 shows the expected result of 0.1x0.1=0.01 sunlit fraction. Also tested with other transmittances of 0.9 and 0.55.

Preliminary fix for overlapping partly transparent shading

c0227ed

mjwitte added the Defect Includes code to repair a defect in EnergyPlus label May 23, 2023

Revised fix for overlapping partly transparent shading

71430dc

mjwitte changed the title ~~Preliminary fix for overlapping partly transparent shading~~ Fix shadow calculations for shading surfaces with transmittances >0 and <1 May 26, 2023

mjwitte changed the title ~~Fix shadow calculations for shading surfaces with transmittances >0 and <1~~ Fix overlapping shadow calculations for shading surfaces with transmittances >0 and <1 May 26, 2023

mjwitte changed the title ~~Fix overlapping shadow calculations for shading surfaces with transmittances >0 and <1~~ Fix overlapping shadow calculations for shading surface transmittance >0 and <1 May 26, 2023

mjwitte added 3 commits June 2, 2023 12:06

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

5cf1143

…ngSched

Skip some setup if no daylighting points - facilitates solar shading …

6032f3b

…unit tests

SolarShadingTest_PolygonOverlap

0700417

rraustad reviewed Jun 4, 2023

View reviewed changes

mjwitte added 12 commits June 5, 2023 13:55

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

c969219

…ngSched

Fix warnings

d8096fd

Revert daylighting change add call to InitSolarViewFactors

fcd41f3

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

3d7e1bf

…ngSched

More DeterminePolygonOverlap tweaks

9ee09b9

More unit tests SolarShadingTest_PolygonOverlap2 and 3

ef33086

PolygonOverlap unit test repairs

e27757b

Add some missing inits

3f9a49c

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

5eb35dc

…ngSched

Fix warnings

2040e6d

Overlapping shadow doc updates

3e369e4

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

3f4ccbc

…ngSched

mjwitte marked this pull request as ready for review June 7, 2023 16:17

mjwitte commented Jun 7, 2023

View reviewed changes

mjwitte requested review from Myoldmopar, nealkruis, rraustad and vidanovic June 7, 2023 17:08

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

5b83888

…ngSched

rraustad reviewed Jun 7, 2023

View reviewed changes

rraustad reviewed Jun 8, 2023

View reviewed changes

rraustad approved these changes Jun 8, 2023

View reviewed changes

vidanovic approved these changes Jun 8, 2023

View reviewed changes

mjwitte added 4 commits June 8, 2023 11:02

Merge remote-tracking branch 'remotes/origin/develop' into 10022Shadi…

2d39e13

…ngSched

More overlapping shadow doc updates

9d6fe9a

Address doc comments

80b448e

Address doc comments - again

6cd72a3

Myoldmopar merged commit 65f1324 into develop Jun 14, 2023

Myoldmopar deleted the 10022ShadingSched branch June 14, 2023 14:59

mjwitte mentioned this pull request Jul 31, 2023

Solar through multiple partially transmitting shading surfaces is not correct at the receiving surface #6420

Closed

3 tasks

shorowit mentioned this pull request Sep 28, 2023

OpenStudio 3.7.0/EnergyPlus 23.2.0 NREL/OpenStudio-HPXML#1476

Merged

7 tasks

		// Check for exceeding array limits.
		if (NS3 > state.dataSolarShading->MaxHCS) {

		if (NV3 < state.dataSolarShading->MaxHCV && NS3 <= state.dataSolarShading->MaxHCS) {
		if (NV3 < state.dataSolarShading->MaxHCV) {

		if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0)
		if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0) {

		\end{equation}

		The actual $\tau$ of overlapping polygons i and j is the product of $\tau$\(_{i}\) and $\tau$\(_{j}\). For example, if a wall is fully shaded by two partially transparent shades, one with a transmittance of 0.8 and one with a transmittance of 0.5, the resulting sunlit fraction would be $0.8 * 0.5 = 0.4$. Because each shadow is first applied with its own transmittance, the transmittance for the overlap in the summation is the difference between the sum and the product of the shadow transmittances:

Fix overlapping shadow calculations for shading surface transmittance >0 and <1 #10040

Fix overlapping shadow calculations for shading surface transmittance >0 and <1 #10040

Conversation

mjwitte commented May 23, 2023 • edited by rraustad Loading

Pull request overview

Diffs

The theory

Pull Request Author

Reviewer

mjwitte commented May 24, 2023

Choose a reason for hiding this comment

rraustad Jun 4, 2023 • edited Loading

Choose a reason for hiding this comment

rraustad Jun 5, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjwitte commented Jun 7, 2023

mjwitte 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

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 Jun 7, 2023 • edited Loading

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 Jun 8, 2023 • edited Loading

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 left a comment

Choose a reason for hiding this comment

vidanovic commented Jun 8, 2023

vidanovic left a comment

Choose a reason for hiding this comment

rraustad commented Jun 9, 2023

rraustad commented Jun 9, 2023

Myoldmopar commented Jun 14, 2023

Myoldmopar commented Jun 14, 2023

ecoeficiente commented Jul 31, 2023

mjwitte commented Jul 31, 2023

mjwitte commented May 23, 2023 •

edited by rraustad

Loading

rraustad Jun 4, 2023 •

edited

Loading

rraustad Jun 5, 2023 •

edited

Loading

rraustad Jun 7, 2023 •

edited

Loading

rraustad Jun 8, 2023 •

edited

Loading