Vectorize pvlib.bifacial.utils._vf_ground_sky_2d across surface_tilt

docs/sphinx/source/whatsnew/v0.9.5.rst

@@ -20,12 +20,14 @@ Enhancements
   :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
 * Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
   support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
-
 * Added an optional ``model`` parameter to
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa`
   to enable use of the hay-davies sky diffuse irradiance model
   instead of the default isotropic model. (:pull:`1668`)
+* Improved execution speed of the infinite sheds model by 10-25% for
+  simulations with time series surface orientation. (:issue:`1680`, :pull:`1682`)
+
 
 Bug fixes
 ~~~~~~~~~
@@ -74,3 +76,4 @@ Contributors
 * Michael Deceglie (:ghuser:`mdeceglie`)
 * Saurabh Aneja (:ghuser:`spaneja`)
 * John Moseley (:ghuser:`johnMoseleyArray`)
+* Areeba Turabi (:ghuser:`aturabi`)

pvlib/bifacial/infinite_sheds.py

@@ -38,17 +38,10 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
 
     Returns
     -------
-    fgnd_sky : float
+    fgnd_sky : numeric
         Integration of view factor over the length between adjacent, interior
-        rows. [unitless]
-    fz : ndarray
-        Fraction of distance from the previous row to the next row. [unitless]
-    fz_sky : ndarray
-        View factors at discrete points between adjacent, interior rows.
-        [unitless]
-
+        rows.  Shape matches that of ``surface_tilt``. [unitless]
     """
-    # TODO: vectorize over surface_tilt
     # Abuse utils._vf_ground_sky_2d by supplying surface_tilt in place
     # of a signed rotation. This is OK because
     # 1) z span the full distance between 2 rows, and
@@ -57,12 +50,10 @@ def _vf_ground_sky_integ(surface_tilt, surface_azimuth, gcr, height,
     # The VFs to the sky will thus be symmetric around z=0.5
     z = np.linspace(0, 1, npoints)
     rotation = np.atleast_1d(surface_tilt)
-    fz_sky = np.zeros((len(rotation), npoints))
-    for k, r in enumerate(rotation):
-        vf, _ = utils._vf_ground_sky_2d(z, r, gcr, pitch, height, max_rows)
-        fz_sky[k, :] = vf
+    fz_sky, _ = utils._vf_ground_sky_2d(z, rotation, gcr, pitch, height,
+                                        max_rows)
     # calculate the integrated view factor for all of the ground between rows
-    return np.trapz(fz_sky, z, axis=1)
+    return np.trapz(fz_sky, z, axis=0)
 
 
 def _poa_ground_shadows(poa_ground, f_gnd_beam, df, vf_gnd_sky):
@@ -469,7 +460,7 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith,
         on the surface that is not reflected away. [unitless]
 
     npoints : int, default 100
-        Number of points used to discretize distance along the ground.
+        Number of discretization points for calculating integrated viewfactors.
 
     Returns
     -------
@@ -701,7 +692,7 @@ def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
         etc. A negative value is a reduction in back irradiance. [unitless]
 
     npoints : int, default 100
-        Number of points used to discretize distance along the ground.
+        Number of discretization points for calculating integrated viewfactors.
 
     Returns
     -------

pvlib/bifacial/utils.py

@@ -5,7 +5,6 @@
 import numpy as np
 from pvlib.tools import sind, cosd, tand
 
-
 def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth):
     """
     Tangent of the angle between the zenith vector and the sun vector
@@ -104,7 +103,7 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
         Position on the ground between two rows, as a fraction of the pitch.
         x = 0 corresponds to the point on the ground directly below the
         center point of a row. Positive x is towards the right. [unitless]
-    rotation : float
+    rotation : numeric
         Rotation angle of the row's right edge relative to row center.
         [degree]
     gcr : float
@@ -120,30 +119,36 @@ def _vf_ground_sky_2d(x, rotation, gcr, pitch, height, max_rows=10):
 
     Returns
     -------
-    vf : numeric
-        Fraction of sky dome visible from each point on the ground. [unitless]
+    vf : array
+        Fraction of sky dome visible from each point on the ground.
+        Shape is (len(x), len(rotation)). [unitless]
     wedge_angles : array
         Angles defining each wedge of sky that is blocked by a row. Shape is
-        (2, len(x), 2*max_rows+1). ``wedge_angles[0,:,:]`` is the
-        starting angle of each wedge, ``wedge_angles[1,:,:]`` is the end angle.
-        [degree]
+        (2, len(x), len(rotation), 2*max_rows+1). ``wedge_angles[0,:,:,:]``
+        is the starting angle of each wedge, ``wedge_angles[1,:,:,:]`` is the
+        end angle. [degree]
     """
-    x = np.atleast_1d(x)  # handle float
+    # handle floats:
+    x = np.atleast_1d(x)[:, np.newaxis, np.newaxis]
+    rotation = np.atleast_1d(rotation)[np.newaxis, :, np.newaxis]
     all_k = np.arange(-max_rows, max_rows + 1)
     width = gcr * pitch / 2.
+    distance_to_row_centers = (all_k - x) * pitch
+    dy = width * sind(rotation)
+    dx = width * cosd(rotation)
     # angles from x to right edge of each row
-    a1 = height + width * sind(rotation)
-    b1 = (all_k - x[:, np.newaxis]) * pitch + width * cosd(rotation)
-    phi_1 = np.degrees(np.arctan2(a1, b1))
+    a1 = height + dy
+    b1 = distance_to_row_centers + dx
+    phi_1 = np.arctan2(a1, b1)  # dimensions: (x, rotation, row)
     # angles from x to left edge of each row
-    a2 = height - width * sind(rotation)
-    b2 = (all_k - x[:, np.newaxis]) * pitch - width * cosd(rotation)
-    phi_2 = np.degrees(np.arctan2(a2, b2))
-    phi = np.stack([phi_1, phi_2])
-    swap = phi[0, :, :] > phi[1, :, :]
-    # swap where phi_1 > phi_2 so that phi_1[0,:,:] is the lesser angle
+    a2 = height - dy
+    b2 = distance_to_row_centers - dx
+    phi_2 = np.arctan2(a2, b2)  # dimensions: (x, rotation, row)
+    phi = np.stack([phi_1, phi_2])  # dimensions: (row edge, x, rotation, row)
+    swap = phi_1 > phi_2
+    # swap where phi_1 > phi_2 so that phi[0,:,:,:] is the lesser angle
     phi = np.where(swap, phi[::-1], phi)
     # right edge of next row - left edge of previous row
-    wedge_vfs = 0.5 * (cosd(phi[1, :, 1:]) - cosd(phi[0, :, :-1]))
-    vf = np.sum(np.where(wedge_vfs > 0, wedge_vfs, 0.), axis=1)
-    return vf, phi
+    wedge_vfs = 0.5 * (np.cos(phi[1, :, :, 1:]) - np.cos(phi[0, :, :, :-1]))
+    vf = np.sum(np.clip(wedge_vfs, a_min=0., a_max=None), axis=-1)
+    return vf, np.degrees(phi)

pvlib/tests/bifacial/test_utils.py

@@ -35,7 +35,7 @@ def test_system_fixed_tilt():
     c22 = (-2 - sqr3) / np.sqrt(1.25**2 + (2 + sqr3)**2)  # right edge row 0
     c23 = (0 - sqr3) / np.sqrt(1.25**2 + (0 - sqr3)**2)  # right edge row 1
     vf_2 = 0.5 * (c23 - c22 + c21 - c20)  # vf at point 1
-    vfs_ground_sky = np.array([vf_0, vf_1, vf_2])
+    vfs_ground_sky = np.array([[vf_0], [vf_1], [vf_2]])
     return syst, pts, vfs_ground_sky