Correct infinite_sheds view factor from row to sky and ground; expose vf functions in bifacial.utils

[mikofski]

• Closes is vf_row_sky correct? #1665
• I am familiar with the contributing guidelines
• Tests added
• Updates entries in docs/sphinx/source/reference for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

This corrects what I believe are typos in the view factors from the module surface to the sky and ground.

[mikofski]

Add a gist to show that this is the correct formulation: https://gist.github.com/mikofski/92724afdeaec78ba816bf76a174d8183

[cwhanse]

@mikofski I push to this PR to correct the three failing tests, as I think the analytic solution has the same conceptual error.

[mikofski]

okay with me. I haven't corrected the tests yet until we were thumbs up to move forward. Thanks!

[kandersolar]

In case anyone else is curious how much this changes simulation results, here's a comparison:

Source

import pvlib
import pvfactors
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

times = pd.date_range(f'2020-01-01 00:00', f'2020-12-31 23:59', freq='5min', tz='Etc/GMT+5')
location = pvlib.location.Location(40, -80)
solpos = location.get_solarposition(times)
solpos = solpos.loc[solpos['apparent_elevation'] > 0, :]
times = solpos.index
irrad = location.get_clearsky(times)
dniet = pvlib.irradiance.get_extra_radiation(times)

sat = pvlib.tracking.singleaxis(solpos.apparent_zenith, solpos.azimuth, gcr=0.5, backtrack=True)

infsheds = pvlib.bifacial.infinite_sheds.get_irradiance(
    surface_tilt=sat.surface_tilt, surface_azimuth=sat.surface_azimuth,
    solar_zenith=solpos.apparent_zenith, solar_azimuth=solpos.azimuth,
    gcr=0.5, height=1.5, pitch=4.0,
    ghi=irrad.ghi, dhi=irrad.dhi, dni=irrad.dni, albedo=0.2,
)

# %%

# run the other version in another console:
infsheds2 = pd.read_clipboard()

# %%

fig, axes = plt.subplots(3, 5)
axes = np.ravel(axes)

for column, ax in zip(infsheds.columns, axes):
    x = infsheds[column].values
    y = infsheds2[column].values
    
    mbe = np.mean(y - x)
    
    ax.scatter(x, y, s=1)
    ax.text(0.1, 0.8, f'MBE={mbe:0.01f} W/m$^2$', transform=ax.transAxes)
    
    ax.axline((0, 0), slope=1, c='k')
    ax.set_xlabel('v0.9.4')
    ax.set_ylabel('ef346dd6')
    ax.set_title(column)

[cwhanse]

Suggestions on how to archive a document document describing the calculation of the known result? I can interpret the code in the unit test but that's because I did that calculation in the not-distant past.

[kandersolar]

I put notes like that in gists sometimes, but I don't know what promises github makes about their longevity. You could have archive.org take a snapshot of the gist for a more permanent record. It doesn't seem to capture the notebook content for ipynb gists, but markdown seems to work fine: https://web.archive.org/web/20230217141453/https://gist.github.com/kanderso-nrel/655f996b08c18109effbdad211a8579e

[mikofski]

Would this gist of vfsky calculated numerically using trapezoid rule and S&H's A-10 be useful as an analog test? To at least it is the purest most easily defensible method to calculate the VF, just wish I was smart enuf to solve it symbolically :/

I realize that A-10 (and B-1) is using a flat plane as proxy for sky, but as in B-71, I think the shape of the virtual surface representing the skydome is arbitrary. EG: one could draw an imaginary plane anywhere between the ground and the sky and so long as it extended infinitely to the horizon it would exchange the same radiation (by which I mean the extrinsic total power and not the intrinsic flux or power per unit area which would differ between virtual surfaces and the skydome b/c areas differ).

[mikofski]

Perhaps out of scope for this PR but a clever way to test isotropic diffuse sky on ground is using superposition, because VF’s of all radiation leaving sky must sum to one therefore isotopic minus combined total of diffuse sky incident on both sides of panel must be equal to the incident sky diffuse on the ground. Note this only accounts for the diffuse fraction of GHI or the second part of the expression derived here:

diffuse sky incident on ground = GHI * df*vf_gnd_sky = isotropic - (diffuse sky on front + diffuse sky on back)

so it’s a simple way to check vf_gnd_sky

[kandersolar]

@mikofski I push to this PR to correct the three failing tests, as I think the analytic solution has the same conceptual error.

@cwhanse are updated tests from you still the next step for this PR? Just making sure this PR isn't being accidentally ignored by reviewers.

[cwhanse]

@cwhanse are updated tests from you still the next step for this PR? Just making sure this PR isn't being accidentally ignored by reviewers.

Yes, I am working on them now.

[cwhanse]

@mikofski I believe this line (and similar for unshaded) are in error:

vf_shade_sky_integ = np.trapz(y, x, axis=0)

If the intent is to compute the average view factor, then the integral needs to be multiplied by 1 / f_x (unshaded by 1/(1 - f_x))

Do you agree?

[mikofski]

Do you agree?

Yes! You are 💯 correct! Nice catch! I would be happy to add a test for this. A simple sanity check is to compare against the original linear average assumption.

However if integrating across the entire module length, then no divisor is necessary because it's normalized length is unity. Maybe that was the oversight/confusion?

[mikofski]

The shade/no-shade integral happens 4 places

in _vf_row_sky_integ():

for shaded:

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 148 in 7a2ec9b

vf_shade_sky_integ = np.trapz(y, x, axis=0)

and unshaded:

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 153 in 7a2ec9b

vf_noshade_sky_integ = np.trapz(y, x, axis=0)

in vf_row_ground_integ():

for shaded:

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 302 in 7a2ec9b

vf_shade_ground_integ = np.trapz(y, x, axis=0)

and unshaded:

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 310 in 7a2ec9b

vf_noshade_ground_integ = np.trapz(y, x, axis=0)

question?

However, I question whether these separate shade & unshaded sections of the module are still needed. They are recombined in _poa_sky_diffuse_pv():

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 180 in 7a2ec9b

return dhi * (f_x * vf_shade_sky_integ + (1 - f_x) * vf_noshade_sky_integ)

As you can be seen, the denominator cancels out and one is left with the sum which would be the same as integrating across the entire length of the module from zero to 1 so that the divisor is unity.

Ditto for _poa_ground_pv():

pvlib-python/pvlib/bifacial/infinite_sheds.py

Line 338 in 7a2ec9b

return poa_ground * (f_x * f_gnd_pv_shade + (1 - f_x) * f_gnd_pv_noshade)

(Note: I believe the docstrings in _poa_ground_pv() should be updated to state "fraction of irradiance from ground incident on shaded|unshaded part of PV surface" as this is how it is applied in get_irradiance_poa())

thoughts?

I'm fine with leaving it as is, but I thought a note would be useful for explaining why it is this way. Originally I was not using trapz() to integrate the view factors along the full length of the module surface, but approximating it as an average assuming that it was mostly linear. To reduce errors from non-linearity, the view factor was approximated by splitting the module into the 2 sections: shaded and unshaded. However, if using trapz() then this 2 section approach is probably not needed, and could be removed. I'll leave it to @pvlib/pvlib-core to decide.

[cwhanse]

I can see the math. There's no reason to do more calculations than needed, but I need to wrap my head around that approach, in the context of the beam diffuse irradiance contributions.

Turns out that the integral of the view factor along the slant surface (once correctly formulated) has a closed form solution. So even if we keep the shade/noshade approach we can replace trapz with an exact expression.

I think the function calculating the view factor from a point on the module surface to the sky has utility outside of infinite_sheds, and I'm inclined to move it to bifacial.utils.

[kandersolar]

@cwhanse do you think an analogous analytical solution exists for the averaged view factor from the ground to a single sky wedge? It would be great if we could get rid of the npoints dimension in the ground-sky VF calculation arrays.

[mikofski]

@cwhanse & @kandersolar how can I help push this along?

[cwhanse]

@mikofski your review would be much appreciated. I lost track of who had the next action (I do) on @kandersolar's comments.

[mikofski]

Hi @cwhanse , can you merge from upstream/main and resolve the conflicts? Then I'll take a look. Thanks!

[kandersolar]

A final bit of hand-wringing, otherwise I think LGTM.

Besides the original bug fix, I like having these vf functions exposed, and infinite_sheds.py feels much more approachable now.

[cwhanse]

@kandersolar I think I've addressed your comments.

@mikofski I don't know that my review is meaningful, since I am a co-author. Another set of eyes on this would be great.

[kandersolar]

@kandersolar I think I've addressed your comments.

I think my last review still has two outstanding comments. In case they aren't showing up, here are links: #1666 (comment) and #1666 (comment)

One other thought: perhaps x0, x1 should have default values 0, 1 in the _integ functions?

[cwhanse]

I'm going to review so the GH reminder is cleared.

@kandersolar do you want me to move the whatsnew content to v0.10.0?

CI failure is codecov.

[kandersolar]

We'll have to merge the 0.9.6 and 0.10.0 whatsnew files in the release finalization PR, so either one is fine in the meantime.

[kandersolar]

I'm going to go ahead and merge this for v0.10.0 since Cliff and I think it is correct. @mikofski don't let that stop you from taking a look too, of course.

[mikofski]

@kevinsa5 and @cwhanse I did a simple test, using the original formulas instead of _vf_poly and I get identical output:

test vf_row_sky_2d():

import pvlib as pvl
import numpy as np
from matplotlib import pyplot as plt

# 100 points on the module surface from bottom (0) to top (1)
x = np.linspace(0, 1, 100)
psix = pvl.shading.masking_angle(surface_tilt, gcr, x)

# original formula according to reference
vfx = 0.5 * (pvl.tools.cosd(psix + surface_tilt) + 1)

# test using new function
vf_test = pvl.bifacial.utils.vf_row_sky_2d(surface_tilt, gcr, x)

np.allclose(vf_test, vfx)
# True

plt.ion()
plt.plot(vfx, vf_test)

test vf_row_ground_2d():

def _ground_angle(x, surface_tilt, gcr):
    """
    Angle from horizontal of the line from a point x on the row slant length
    to the bottom of the facing row.

    The angles are clockwise from horizontal, rather than the usual
    counterclockwise direction.

    Parameters
    ----------
    x : numeric
        fraction of row slant length from bottom, ``x = 0`` is at the row
        bottom, ``x = 1`` is at the top of the row.
    surface_tilt : numeric
        Surface tilt angle in degrees from horizontal, e.g., surface facing up
        = 0, surface facing horizon = 90. [degree]
    gcr : float
        ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]

    Returns
    -------
    psi : numeric
        Angle [degree].
    """
    #  : \\            \
    #  :  \\            \
    #  :   \\            \
    #  :    \\            \  facing row
    #  :     \\.___________\
    #  :       \  ^*-.  psi \
    #  :        \  x   *-.   \
    #  :         \  v      *-.\
    #  :          \<-----P---->\

    x1 = gcr * x * pvl.tools.sind(surface_tilt)
    x2 = gcr * x * pvl.tools.cosd(surface_tilt) + 1
    psi = np.arctan2(x1, x2)  # do this first because it handles 0 / 0
    return np.rad2deg(psi)


psib = _ground_angle(x, surface_tilt, gcr)

# original formaul from reference
vf_gnd_x = 0.5 * (1 - pvl.tools.cosd(psib - surface_tilt))

# new formula
vf_gnd_test = pvl.bifacial.utils.vf_row_ground_2d(surface_tilt, gcr, x)

plt.plot(vf_gnd_x, vf_gnd_test)
plt.xlabel('original formula from reference')
plt.ylabel('new formulation')
plt.title('VF from point "x" on module surface to ground')

Personally I find the new formulas slightly opaque and I miss the LaTeX:

original formula for vf_row_sky() from reference

$$\mathit{VF}_{row \rightarrow sky} = \frac{1 + \cos\left(\psi_t + \beta\right)}{2}$$

original formula for vf_row_gnd() from reference

$$\mathit{VF}_{row \rightarrow ground} = \frac{1 + \cos\left(\psi_b - \beta\right)}{2}$$