From 162027e2213ae2a6cdcee552eef95df9a7179527 Mon Sep 17 00:00:00 2001 From: "Marc A. Anoma" Date: Wed, 20 Nov 2019 16:34:39 +0900 Subject: [PATCH] Merge new AOI methods into full mode workflow (#94) Fixes #36 * Add workflow for removing aoi losses and calculating q_absorbed * The engine needs to fit vf calculator as well * Speed up calculation by ignoring "null" surfaces * ~2x speed up * Add "is_empty" property to surfaces * Speed up vf methods by ignoring null surfaces * Use "is_empty" flag in aoi methods as well * Test engine calculated absorbed values with diffuse AOI losses * Check that VF calculator uses reflectivity matrix for diffuse faoi * Test for pv engine using faoi functions * Update both run functions to accept vf calculator params * run_timeseries_engine works but problem with run_timeseries_parallel * Fix so that can run faoi_functions with parallel mode * Try to fix error in Python 2 --- pvfactors/engine.py | 28 +++- pvfactors/geometry/pvarray.py | 1 + pvfactors/geometry/timeseries.py | 6 + pvfactors/irradiance/base.py | 2 +- pvfactors/run.py | 32 ++-- pvfactors/tests/conftest.py | 1 + pvfactors/tests/test_engine.py | 63 +++++++- pvfactors/tests/test_run.py | 148 ++++++++++++++++++ .../tests/test_viewfactors/test_calculator.py | 32 +++- pvfactors/viewfactors/aoimethods.py | 126 ++++++++------- pvfactors/viewfactors/calculator.py | 59 +++++-- pvfactors/viewfactors/vfmethods.py | 30 ++++ 12 files changed, 435 insertions(+), 93 deletions(-) diff --git a/pvfactors/engine.py b/pvfactors/engine.py index 975303a7..41c438a4 100644 --- a/pvfactors/engine.py +++ b/pvfactors/engine.py @@ -164,6 +164,9 @@ def fit(self, timestamps, DNI, DHI, solar_zenith, solar_azimuth, # Add timeseries irradiance results to pvarray self.irradiance.transform(self.pvarray) + # Fit VF calculator + self.vf_calculator.fit(self.n_points) + # Skip timesteps when: # - solar zenith > 90, ie the sun is down # - DNI or DHI is negative, which does not make sense @@ -193,10 +196,10 @@ def run_full_mode(self, fn_build_report=None): # Get the irradiance modeling matrices # shape = n_surfaces, n_timesteps - irradiance_mat, _, invrho_mat, _ = \ + irradiance_mat, rho_mat, invrho_mat, _ = \ self.irradiance.get_full_ts_modeling_vectors(pvarray) - # Calculate view factors + # --- Calculate view factors # shape = n_surfaces, n_surfaces, n_timesteps ts_vf_matrix = self.vf_calculator.build_ts_vf_matrix(pvarray) pvarray.ts_vf_matrix = ts_vf_matrix @@ -204,6 +207,7 @@ def run_full_mode(self, fn_build_report=None): # shape = n_timesteps, n_surfaces, n_surfaces ts_vf_matrix_reshaped = np.moveaxis(ts_vf_matrix, -1, 0) + # --- Solve mathematical problem # Build matrix of inverse reflectivities # shape = n_surfaces, n_surfaces invrho_mat = np.diag(invrho_mat[:, 0]) @@ -220,18 +224,32 @@ def run_full_mode(self, fn_build_report=None): # shape = n_surfaces, n_timesteps qinc = np.dot(invrho_mat, q0) - # Derive other irradiance terms + # --- Derive other irradiance terms # shape = n_surfaces, n_timesteps isotropic_mat = ts_vf_matrix[:-1, -1, :] * irradiance_mat[-1, :] reflection_mat = qinc[:-1, :] - irradiance_mat[:-1, :] - isotropic_mat - # Update surfaces with values: the list is ordered by index + # --- Calculate AOI losses and absorbed irradiance + rho_mat = np.tile(rho_mat[:, 0], (rho_mat.shape[0], 1)).T + # shape [n_surfaces + 1, n_surfaces + 1, n_timestamps] + vf_aoi_matrix = (self.vf_calculator + .build_ts_vf_aoi_matrix(pvarray, rho_mat)) + pvarray.ts_vf_aoi_matrix = vf_aoi_matrix + # shape [n_surfaces, n_surfaces] + irradiance_abs_mat = ( + self.irradiance.get_summed_components(pvarray, absorbed=True)) + # Calculate absorbed irradiance + qabs = (np.einsum('ijk,jk->ik', vf_aoi_matrix, q0)[:-1, :] + + irradiance_abs_mat) + + # --- Update surfaces with values: the lists are ordered by index for idx_surf, ts_surface in enumerate(pvarray.all_ts_surfaces): ts_surface.update_params( {'q0': q0[idx_surf, :], 'qinc': qinc[idx_surf, :], 'isotropic': isotropic_mat[idx_surf, :], - 'reflection': reflection_mat[idx_surf, :]}) + 'reflection': reflection_mat[idx_surf, :], + 'qabs': qabs[idx_surf, :]}) # Return report if function was passed report = None if fn_build_report is None else fn_build_report(pvarray) diff --git a/pvfactors/geometry/pvarray.py b/pvfactors/geometry/pvarray.py index dae33659..a8cecbba 100644 --- a/pvfactors/geometry/pvarray.py +++ b/pvfactors/geometry/pvarray.py @@ -69,6 +69,7 @@ def __init__(self, axis_azimuth=None, gcr=None, pvrow_height=None, # These attributes will be updated by the engine self.ts_vf_matrix = None + self.ts_vf_aoi_matrix = None @classmethod def init_from_dict(cls, pvarray_params, param_names=None): diff --git a/pvfactors/geometry/timeseries.py b/pvfactors/geometry/timeseries.py index 70d83d41..863f147b 100644 --- a/pvfactors/geometry/timeseries.py +++ b/pvfactors/geometry/timeseries.py @@ -249,6 +249,12 @@ def u_vector(self): np.array([-self.n_vector[1, :], self.n_vector[0, :]])) return u_vector + @property + def is_empty(self): + """Check if surface is "empty" by checking if its length is always + zero""" + return np.nansum(self.length) < DISTANCE_TOLERANCE + class TsLineCoords(object): """Timeseries line coordinates class: will provide a helpful shapely-like diff --git a/pvfactors/irradiance/base.py b/pvfactors/irradiance/base.py index b26d3167..00d37867 100644 --- a/pvfactors/irradiance/base.py +++ b/pvfactors/irradiance/base.py @@ -118,7 +118,7 @@ def get_summed_components(self, pvarray, absorbed=True): for component in list_components: value += ts_surface.get_param(component) irradiance_mat.append(value) - return irradiance_mat + return np.array(irradiance_mat) def update_ts_surface_sky_term(self, ts_surface, name_sky_term='sky_term'): """Update the 'sky_term' parameter of a timeseries surface. diff --git a/pvfactors/run.py b/pvfactors/run.py index 69563618..cc14de27 100644 --- a/pvfactors/run.py +++ b/pvfactors/run.py @@ -25,6 +25,7 @@ def run_timeseries_engine(fn_build_report, pvarray_parameters, fast_mode_pvrow_index=None, fast_mode_segment_index=None, irradiance_model_params=None, + vf_calculator_params=None, ghi=None): """Run timeseries simulation without multiprocessing. This is the functional approach to the :py:class:`~pvfactors.engine.PVEngine` class. @@ -77,6 +78,9 @@ def run_timeseries_engine(fn_build_report, pvarray_parameters, irradiance_model_params : dict, optional Dictionary of parameters that will be passed to the irradiance model class as kwargs at instantiation (Default = None) + vf_calculator_params : dict, optional + Dictionary of parameters that will be passed to the VF calculator + class as kwargs at instantiation (Default = None) ghi : array-like, optional Global horizontal irradiance values [W/m2] (Default = None) @@ -87,18 +91,17 @@ class as kwargs at instantiation (Default = None) function """ - # Prepare irradiance model inputs - irradiance_model_params = ({} if irradiance_model_params is None - else irradiance_model_params) + # Prepare input parameters + irradiance_model_params = irradiance_model_params or {} + vf_calculator_params = vf_calculator_params or {} # Instantiate classes and engine irradiance_model = cls_irradiance(**irradiance_model_params) + vf_calculator = cls_vf(**vf_calculator_params) pvarray = cls_pvarray.init_from_dict(pvarray_parameters) - vf_calculator = cls_vf() eng = cls_engine(pvarray, irradiance_model=irradiance_model, vf_calculator=vf_calculator, fast_mode_pvrow_index=fast_mode_pvrow_index, fast_mode_segment_index=fast_mode_segment_index) - # Fit engine eng.fit(timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, surface_azimuth, albedo, ghi=ghi) @@ -118,7 +121,8 @@ def run_parallel_engine(report_builder, pvarray_parameters, cls_irradiance=HybridPerezOrdered, cls_vf=VFCalculator, fast_mode_pvrow_index=None, fast_mode_segment_index=None, - irradiance_model_params=None, n_processes=2, + irradiance_model_params=None, + vf_calculator_params=None, n_processes=2, ghi=None): """Run timeseries simulation using multiprocessing. Here, instead of a function that will build the report, the users will need to pass a class @@ -173,6 +177,9 @@ def run_parallel_engine(report_builder, pvarray_parameters, irradiance_model_params : dict, optional Dictionary of parameters that will be passed to the irradiance model class as kwargs at instantiation (Default = None) + vf_calculator_params : dict, optional + Dictionary of parameters that will be passed to the VF calculator + class as kwargs at instantiation (Default = None) n_processes : int, optional Number of parallel processes to run for the calculation (Default = 2) ghi : array-like, optional @@ -192,9 +199,9 @@ class (or object) if np.isscalar(albedo): albedo = albedo * np.ones(len(dni)) - # Prepare irradiance model inputs - irradiance_model_params = ({} if irradiance_model_params is None - else irradiance_model_params) + # Prepare class input parameters + irradiance_model_params = irradiance_model_params or {} + vf_calculator_params = vf_calculator_params or {} # Fix: np.array_split doesn't work well on pd.DatetimeIndex objects if isinstance(timestamps, pd.DatetimeIndex): @@ -216,6 +223,7 @@ class (or object) folds_fast_mode_pvrow_index = [fast_mode_pvrow_index] * n_processes folds_fast_mode_segment_index = [fast_mode_segment_index] * n_processes folds_irradiance_model_params = [irradiance_model_params] * n_processes + folds_vf_calculator_params = [vf_calculator_params] * n_processes folds_ghi = ([ghi] * n_processes if ghi is None else np.array_split(ghi, n_processes)) report_indices = list(range(n_processes)) @@ -228,7 +236,8 @@ class (or object) folds_cls_engine, folds_cls_irradiance, folds_cls_vf, folds_fast_mode_pvrow_index, folds_fast_mode_segment_index, - folds_irradiance_model_params, folds_ghi, + folds_irradiance_model_params, + folds_vf_calculator_params, folds_ghi, report_indices)) # Start multiprocessing @@ -272,7 +281,7 @@ class (or object) solar_zenith, solar_azimuth, surface_tilt, surface_azimuth,\ albedo, cls_pvarray, cls_engine, cls_irradiance, cls_vf, \ fast_mode_pvrow_index, fast_mode_segment_index, \ - irradiance_model_params, ghi, idx = args + irradiance_model_params, vf_calculator_params, ghi, idx = args report = run_timeseries_engine( report_builder.build, pvarray_parameters, @@ -283,6 +292,7 @@ class (or object) fast_mode_pvrow_index=fast_mode_pvrow_index, fast_mode_segment_index=fast_mode_segment_index, irradiance_model_params=irradiance_model_params, + vf_calculator_params=vf_calculator_params, ghi=ghi) return report, idx diff --git a/pvfactors/tests/conftest.py b/pvfactors/tests/conftest.py index a7508f54..e4028544 100644 --- a/pvfactors/tests/conftest.py +++ b/pvfactors/tests/conftest.py @@ -138,6 +138,7 @@ def fn_report_example(): def fn_report(pvarray): return { 'qinc_front': pvarray.ts_pvrows[1].front.get_param_weighted('qinc'), 'qinc_back': pvarray.ts_pvrows[1].back.get_param_weighted('qinc'), + 'qabs_back': pvarray.ts_pvrows[1].back.get_param_weighted('qabs'), 'iso_front': pvarray.ts_pvrows[1] .front.get_param_weighted('isotropic'), 'iso_back': pvarray.ts_pvrows[1].back.get_param_weighted('isotropic')} diff --git a/pvfactors/tests/test_engine.py b/pvfactors/tests/test_engine.py index 809d720e..254b2b5d 100644 --- a/pvfactors/tests/test_engine.py +++ b/pvfactors/tests/test_engine.py @@ -41,6 +41,10 @@ def test_pvengine_float_inputs_iso(params): pvarray.ts_pvrows[1].front.get_param_weighted('qinc'), 1099.6948573) np.testing.assert_almost_equal( pvarray.ts_pvrows[2].front.get_param_weighted('qinc'), 1102.76149246) + # Check absorbed + np.testing.assert_almost_equal( + pvarray.ts_pvrows[1].front.get_param_weighted('qabs'), + 1099.6948573 * 0.99) def test_pvengine_float_inputs_perez(params): @@ -80,6 +84,13 @@ def test_pvengine_float_inputs_perez(params): np.testing.assert_almost_equal( pvarray.ts_pvrows[1].back.get_param_weighted('qinc'), 116.49050349491208) + # Check absorbed irradiance + np.testing.assert_almost_equal( + pvarray.ts_pvrows[2].front.get_param_weighted('qabs'), + 1112.37717553 * 0.99) + np.testing.assert_almost_equal( + pvarray.ts_pvrows[1].back.get_param_weighted('qabs'), + 116.49050349491208 * 0.97) def test_pvengine_ts_inputs_perez(params_serial, @@ -115,6 +126,8 @@ def test_pvengine_ts_inputs_perez(params_serial, report['iso_front'], [42.816637, 42.780206]) np.testing.assert_array_almost_equal( report['iso_back'], [1.727308, 1.726535]) + np.testing.assert_array_almost_equal( + report['qabs_back'], report['qinc_back'] * 0.97) def test_run_fast_mode_isotropic(params): @@ -545,7 +558,7 @@ def test_check_direct_shading_continuity(): np.testing.assert_allclose(out, expected_out) -def test_create_with_rho_init(params, pvmodule_canadian): +def test_create_engine_with_rho_init(params, pvmodule_canadian): """Check that can create PV engine with rho initialization from faoi functions""" # Create inputs @@ -559,3 +572,51 @@ def test_create_with_rho_init(params, pvmodule_canadian): # Check that rho values are the ones calculated np.testing.assert_allclose(engine.irradiance.rho_front, 0.02900688) np.testing.assert_allclose(engine.irradiance.rho_back, 0.02900688) + + +def test_create_engine_with_rho_init(params, pvmodule_canadian): + """Run PV engine calcs with faoi functions for AOI losses""" + + # Irradiance inputs + timestamps = dt.datetime(2019, 6, 11, 11) + DNI = 1000. + DHI = 100. + + irradiance_model = HybridPerezOrdered() + pvarray = OrderedPVArray.init_from_dict(params) + faoi_fn = faoi_fn_from_pvlib_sandia(pvmodule_canadian) + vfcalculator = VFCalculator(faoi_fn_front=faoi_fn, faoi_fn_back=faoi_fn) + eng = PVEngine(pvarray, irradiance_model=irradiance_model, + vf_calculator=vfcalculator) + + # Make sure aoi methods are available + assert eng.vf_calculator.vf_aoi_methods is not None + + # Fit engine + eng.fit(timestamps, DNI, DHI, + params['solar_zenith'], + params['solar_azimuth'], + params['surface_tilt'], + params['surface_azimuth'], + params['rho_ground']) + + # Run timestep + pvarray = eng.run_full_mode(fn_build_report=lambda pvarray: pvarray) + # Checks + np.testing.assert_almost_equal( + pvarray.ts_pvrows[0].front.get_param_weighted('qinc'), + 1110.1164773159298) + np.testing.assert_almost_equal( + pvarray.ts_pvrows[1].front.get_param_weighted('qinc'), 1110.595903991) + np.testing.assert_almost_equal( + pvarray.ts_pvrows[2].front.get_param_weighted('qinc'), 1112.37717553) + np.testing.assert_almost_equal( + pvarray.ts_pvrows[1].back.get_param_weighted('qinc'), + 116.49050349491208) + # Check absorbed irradiance: calculated using faoi functions + np.testing.assert_almost_equal( + pvarray.ts_pvrows[2].front.get_param_weighted('qabs'), + [1099.1251094]) + np.testing.assert_almost_equal( + pvarray.ts_pvrows[1].back.get_param_weighted('qabs'), + [114.4690984]) diff --git a/pvfactors/tests/test_run.py b/pvfactors/tests/test_run.py index 28370e59..25c10eaa 100644 --- a/pvfactors/tests/test_run.py +++ b/pvfactors/tests/test_run.py @@ -1,5 +1,7 @@ from pvfactors.run import run_timeseries_engine, run_parallel_engine from pvfactors.report import ExampleReportBuilder +from pvfactors.viewfactors.aoimethods import faoi_fn_from_pvlib_sandia +from pvfactors.viewfactors.calculator import VFCalculator import numpy as np import mock @@ -211,3 +213,149 @@ def merge(reports): report[key] += other_report[key] return report + + +def test_run_timeseries_faoi_fn(params_serial, pvmodule_canadian, + df_inputs_clearsky_8760): + """Test that in run_timeseries function, faoi functions are used + correctly""" + # Prepare timeseries inputs + df_inputs = df_inputs_clearsky_8760.iloc[:24, :] + timestamps = df_inputs.index + dni = df_inputs.dni.values + dhi = df_inputs.dhi.values + solar_zenith = df_inputs.solar_zenith.values + solar_azimuth = df_inputs.solar_azimuth.values + surface_tilt = df_inputs.surface_tilt.values + surface_azimuth = df_inputs.surface_azimuth.values + + expected_qinc = 542.018551 + + # --- Test without passing vf parameters + # report function with test in it + def report_fn_with_tests_no_faoi(pvarray): + vf_aoi_matrix = pvarray.ts_vf_aoi_matrix + pvrow = pvarray.ts_pvrows[0] + list_back_pvrow_idx = [ts_surf.index for ts_surf + in pvarray.ts_pvrows[0].all_ts_surfaces] + # Check that sum of vf_aoi is equal to reflectivity values + # since no faoi_fn used + np.testing.assert_allclose( + vf_aoi_matrix[list_back_pvrow_idx, :, 12].sum(axis=1), + [0.99, 0., 0.97, 0.]) + + return {'qinc_back': pvrow.back.get_param_weighted('qinc'), + 'qabs_back': pvrow.back.get_param_weighted('qabs')} + + # create calculator + report = run_timeseries_engine( + report_fn_with_tests_no_faoi, params_serial, + timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, + surface_azimuth, params_serial['rho_ground'], + vf_calculator_params=None) + + np.testing.assert_allclose(np.nansum(report['qinc_back']), expected_qinc) + np.testing.assert_allclose(np.nansum(report['qabs_back']), 525.757995) + + # --- Test when passing vf parameters + # Prepare vf calc params + faoi_fn = faoi_fn_from_pvlib_sandia(pvmodule_canadian) + # the following is a very high number to get agreement in + # integral sums between back and front surfaces + n_sections = 10000 + vf_calc_params = {'faoi_fn_front': faoi_fn, + 'faoi_fn_back': faoi_fn, + 'n_aoi_integral_sections': n_sections} + + def report_fn_with_tests_w_faoi(pvarray): + vf_aoi_matrix = pvarray.ts_vf_aoi_matrix + pvrow = pvarray.ts_pvrows[0] + + list_back_pvrow_idx = [ts_surf.index for ts_surf + in pvrow.all_ts_surfaces] + # Check that sum of vf_aoi is consistent + np.testing.assert_allclose( + vf_aoi_matrix[list_back_pvrow_idx, :, 12].sum(axis=1), + [0.97102, 0., 0.971548, 0.], atol=0, rtol=1e-6) + + return {'qinc_back': pvrow.back.get_param_weighted('qinc'), + 'qabs_back': pvrow.back.get_param_weighted('qabs')} + # create calculator + report = run_timeseries_engine( + report_fn_with_tests_w_faoi, params_serial, + timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, + surface_azimuth, params_serial['rho_ground'], + vf_calculator_params=vf_calc_params) + + np.testing.assert_allclose(np.nansum(report['qinc_back']), expected_qinc) + np.testing.assert_allclose(np.nansum(report['qabs_back']), 522.299276) + + +def test_run_parallel_faoi_fn(params_serial, df_inputs_clearsky_8760): + """Test that in run_parallel function, faoi functions are used + correctly""" + # Prepare timeseries inputs + df_inputs = df_inputs_clearsky_8760.iloc[:24, :] + timestamps = df_inputs.index + dni = df_inputs.dni.values + dhi = df_inputs.dhi.values + solar_zenith = df_inputs.solar_zenith.values + solar_azimuth = df_inputs.solar_azimuth.values + surface_tilt = df_inputs.surface_tilt.values + surface_azimuth = df_inputs.surface_azimuth.values + + expected_qinc = 542.018551 + # create calculator + report = run_parallel_engine( + TestFAOIReportBuilder, params_serial, + timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, + surface_azimuth, params_serial['rho_ground'], + vf_calculator_params=None) + + np.testing.assert_allclose(np.nansum(report['qinc_back']), expected_qinc) + np.testing.assert_allclose(np.nansum(report['qabs_back']), 525.757995) + + # --- Test when passing vf parameters + # the following is a very high number to get agreement in + # integral sums between back and front surfaces + n_sections = 10000 + vf_calc_params = {'faoi_fn_front': FaoiClass, + 'faoi_fn_back': FaoiClass, + 'n_aoi_integral_sections': n_sections} + + # create calculator + report = run_parallel_engine( + TestFAOIReportBuilder, params_serial, + timestamps, dni, dhi, solar_zenith, solar_azimuth, surface_tilt, + surface_azimuth, params_serial['rho_ground'], + vf_calculator_params=vf_calc_params) + + np.testing.assert_allclose(np.nansum(report['qinc_back']), expected_qinc) + np.testing.assert_allclose(np.nansum(report['qabs_back']), 522.299276) + + +class TestFAOIReportBuilder(object): + + @staticmethod + def build(pvarray): + return {'qinc_back': pvarray.ts_pvrows[0].back + .get_param_weighted('qinc').tolist(), + 'qabs_back': pvarray.ts_pvrows[0].back + .get_param_weighted('qabs').tolist()} + + @staticmethod + def merge(reports): + report = reports[0] + keys = report.keys() + for other_report in reports[1:]: + for key in keys: + report[key] += other_report[key] + return report + + +class FaoiClass(object): + + @staticmethod + def faoi(*args, **kwargs): + fn = faoi_fn_from_pvlib_sandia('Canadian_Solar_CS5P_220M___2009_') + return fn(*args, **kwargs) diff --git a/pvfactors/tests/test_viewfactors/test_calculator.py b/pvfactors/tests/test_viewfactors/test_calculator.py index ba3c21a7..324116cf 100644 --- a/pvfactors/tests/test_viewfactors/test_calculator.py +++ b/pvfactors/tests/test_viewfactors/test_calculator.py @@ -33,7 +33,8 @@ def test_vfcalculator_vectorized(params): def test_vfcalculator_aoi_methods(params): - """Check that calculation of vf_aoi_matrix makes sense: + """Check that calculation of vf_aoi_matrix makes sense when using + faoi function is passed: all NREL-like vfs should sum up to 1 when taking many integration points""" # Prepare pv array @@ -49,7 +50,7 @@ def faoi_fn(aoi_angles): return np.ones_like(aoi_angles) vfcalculator = VFCalculator( faoi_fn, faoi_fn, n_aoi_integral_sections=n_integration_sections) vfcalculator.fit(n_timestamps) - vf_aoi_matrix = vfcalculator.build_ts_vf_aoi_matrix(pvarray) + vf_aoi_matrix = vfcalculator.build_ts_vf_aoi_matrix(pvarray, None) # Check that correct size assert vf_aoi_matrix.shape == (47, 47, 1) @@ -68,6 +69,33 @@ def faoi_fn(aoi_angles): return np.ones_like(aoi_angles) atol=0, rtol=1.1e-3) +def test_vfcalculator_no_aoi_functions(params): + """Check that the VF calculator calculates aoi losses + with diffuse reflection values when no faoi_fn is passed""" + + # Prepare pv array + params.update({'cut': {0: {'front': 3}, 1: {'back': 2}}}) + pvarray = OrderedPVArray.fit_from_dict_of_scalars(params) + n_timestamps = 1 + # Create calculator without faoi functions + vfcalculator = VFCalculator() + vfcalculator.fit(n_timestamps) + # Run calculation of view factors first + vf_matrix = vfcalculator.build_ts_vf_matrix(pvarray) + # Make sure that the matrix was saved + assert np.sum(vfcalculator.vf_matrix) > 0 + # Compute reflectivity + rho_mat = np.tile([0.03] * (pvarray.n_ts_surfaces + 1), + (pvarray.n_ts_surfaces + 1, 1)).T + assert rho_mat.shape == (pvarray.n_ts_surfaces + 1, + pvarray.n_ts_surfaces + 1) + vf_aoi_matrix = vfcalculator.build_ts_vf_aoi_matrix(pvarray, rho_mat) + + # Check that correct size + assert vf_aoi_matrix.shape == (47, 47, 1) + np.testing.assert_allclose(vf_matrix * 0.97, vf_aoi_matrix) + + def test_ts_view_factors(): """Test calculation of timeseries view factors for center PV row""" # Create base params diff --git a/pvfactors/viewfactors/aoimethods.py b/pvfactors/viewfactors/aoimethods.py index c847a543..bb460252 100644 --- a/pvfactors/viewfactors/aoimethods.py +++ b/pvfactors/viewfactors/aoimethods.py @@ -294,31 +294,37 @@ def _vf_aoi_surface_to_surface(self, surf_1, surf_2, is_back=True): is_back : bool Flag specifying whether pv row surface is on back or front side of PV row (Default = True) + Returns ------- vf_aoi : np.ndarray View factors with aoi losses from surface 1 to surface 2, dimension is [n_timesteps] """ - # Get surface 1 params - u_vector = surf_1.u_vector - centroid = surf_1.centroid - # Calculate AOI angles - aoi_angles_1 = self._calculate_aoi_angles(u_vector, centroid, - surf_2.b1) - aoi_angles_2 = self._calculate_aoi_angles(u_vector, centroid, - surf_2.b2) - low_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_1, - aoi_angles_2) - high_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_2, - aoi_angles_1) - # Calculate vf_aoi - vf_aoi_raw = self._calculate_vf_aoi_wedge_level( - low_aoi_angles, high_aoi_angles, is_back=is_back) - # Should be zero where either of the surfaces have zero length - vf_aoi = np.where((surf_1.length < DISTANCE_TOLERANCE) - | (surf_2.length < DISTANCE_TOLERANCE), 0., - vf_aoi_raw) + # skip calculation if either surface is empty (always zero length) + skip = surf_1.is_empty or surf_2.is_empty + if skip: + vf_aoi = np.zeros_like(surf_2.length) + else: + # Get surface 1 params + u_vector = surf_1.u_vector + centroid = surf_1.centroid + # Calculate AOI angles + aoi_angles_1 = self._calculate_aoi_angles(u_vector, centroid, + surf_2.b1) + aoi_angles_2 = self._calculate_aoi_angles(u_vector, centroid, + surf_2.b2) + low_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_1, + aoi_angles_2) + high_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_2, + aoi_angles_1) + # Calculate vf_aoi + vf_aoi_raw = self._calculate_vf_aoi_wedge_level( + low_aoi_angles, high_aoi_angles, is_back=is_back) + # Should be zero where either of the surfaces have zero length + vf_aoi = np.where((surf_1.length < DISTANCE_TOLERANCE) + | (surf_2.length < DISTANCE_TOLERANCE), 0., + vf_aoi_raw) return vf_aoi def _vf_aoi_pvrow_surf_to_gnd_surf_obstruction( @@ -363,45 +369,51 @@ def _vf_aoi_pvrow_surf_to_gnd_surf_obstruction( View factors aoi from timeseries PV row surface to timeseries ground surface, dimension is [n_timesteps] """ - centroid = pvrow_surf.centroid - u_vector = pvrow_surf.u_vector - no_obstruction = (is_left & (pvrow_idx == 0)) \ - or ((not is_left) & (pvrow_idx == n_pvrows - 1)) - if no_obstruction: - # There is no obstruction to the ground surface - aoi_angles_1 = self._calculate_aoi_angles(u_vector, centroid, - gnd_surf.b1) - aoi_angles_2 = self._calculate_aoi_angles(u_vector, centroid, - gnd_surf.b2) - else: - # Get lowest point of obstructing point - idx_obstructing_pvrow = pvrow_idx - 1 if is_left else pvrow_idx + 1 - pt_obstr = ts_pvrows[idx_obstructing_pvrow - ].full_pvrow_coords.lowest_point - # Adjust angle seen when there is obstruction - aoi_angles_1 = self._calculate_aoi_angles_w_obstruction( - u_vector, centroid, gnd_surf.b1, pt_obstr, is_left) - aoi_angles_2 = self._calculate_aoi_angles_w_obstruction( - u_vector, centroid, gnd_surf.b2, pt_obstr, is_left) - - low_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_1, - aoi_angles_2) - high_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, aoi_angles_2, - aoi_angles_1) - vf_aoi_raw = self._calculate_vf_aoi_wedge_level( - low_aoi_angles, high_aoi_angles, is_back=is_back) - # Should be zero where either of the surfaces have zero length - vf_aoi_raw = np.where((ts_length < DISTANCE_TOLERANCE) - | (gnd_surf.length < DISTANCE_TOLERANCE), 0., - vf_aoi_raw) - - # Final result depends on whether front or back surface - if is_left: - vf_aoi = (np.where(tilted_to_left, 0., vf_aoi_raw) if is_back - else np.where(tilted_to_left, vf_aoi_raw, 0.)) + # skip calculation if either surface is empty (always zero length) + skip = pvrow_surf.is_empty or gnd_surf.is_empty + if skip: + vf_aoi = np.zeros_like(gnd_surf.length) else: - vf_aoi = (np.where(tilted_to_left, vf_aoi_raw, 0.) if is_back - else np.where(tilted_to_left, 0., vf_aoi_raw)) + centroid = pvrow_surf.centroid + u_vector = pvrow_surf.u_vector + no_obstruction = (is_left & (pvrow_idx == 0)) \ + or ((not is_left) & (pvrow_idx == n_pvrows - 1)) + if no_obstruction: + # There is no obstruction to the ground surface + aoi_angles_1 = self._calculate_aoi_angles(u_vector, centroid, + gnd_surf.b1) + aoi_angles_2 = self._calculate_aoi_angles(u_vector, centroid, + gnd_surf.b2) + else: + # Get lowest point of obstructing point + idx_obstructing_pvrow = (pvrow_idx - 1 if is_left + else pvrow_idx + 1) + pt_obstr = ts_pvrows[idx_obstructing_pvrow + ].full_pvrow_coords.lowest_point + # Adjust angle seen when there is obstruction + aoi_angles_1 = self._calculate_aoi_angles_w_obstruction( + u_vector, centroid, gnd_surf.b1, pt_obstr, is_left) + aoi_angles_2 = self._calculate_aoi_angles_w_obstruction( + u_vector, centroid, gnd_surf.b2, pt_obstr, is_left) + + low_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, + aoi_angles_1, aoi_angles_2) + high_aoi_angles = np.where(aoi_angles_1 < aoi_angles_2, + aoi_angles_2, aoi_angles_1) + vf_aoi_raw = self._calculate_vf_aoi_wedge_level( + low_aoi_angles, high_aoi_angles, is_back=is_back) + # Should be zero where either of the surfaces have zero length + vf_aoi_raw = np.where((ts_length < DISTANCE_TOLERANCE) + | (gnd_surf.length < DISTANCE_TOLERANCE), 0., + vf_aoi_raw) + + # Final result depends on whether front or back surface + if is_left: + vf_aoi = (np.where(tilted_to_left, 0., vf_aoi_raw) if is_back + else np.where(tilted_to_left, vf_aoi_raw, 0.)) + else: + vf_aoi = (np.where(tilted_to_left, vf_aoi_raw, 0.) if is_back + else np.where(tilted_to_left, 0., vf_aoi_raw)) return vf_aoi diff --git a/pvfactors/viewfactors/calculator.py b/pvfactors/viewfactors/calculator.py index 62525dba..dc337d4e 100644 --- a/pvfactors/viewfactors/calculator.py +++ b/pvfactors/viewfactors/calculator.py @@ -20,13 +20,15 @@ def __init__(self, faoi_fn_front=None, faoi_fn_back=None, Parameters ---------- - faoi_fn_front : function, optional - Function which takes a list (or numpy array) of incidence angles + faoi_fn_front : function or object, optional + Function (or object containing ``faoi`` method) + which takes a list (or numpy array) of incidence angles measured from the surface horizontal (with values from 0 to 180 deg) and returns the fAOI values for the front side of PV rows (default = None) - faoi_fn_back : function, optional - Function which takes a list (or numpy array) of incidence angles + faoi_fn_back : function or object, optional + Function (or object containing ``faoi`` method) + which takes a list (or numpy array) of incidence angles measured from the surface horizontal (with values from 0 to 180 deg) and returns the fAOI values for the back side of PV rows (default = None) @@ -39,6 +41,13 @@ def __init__(self, faoi_fn_front=None, faoi_fn_back=None, if (faoi_fn_front is None) or (faoi_fn_back is None): self.vf_aoi_methods = None else: + # Check whether got function or object, and take ``faoi`` method + # if object was passed + faoi_fn_front = (faoi_fn_front.faoi + if hasattr(faoi_fn_front, 'faoi') + else faoi_fn_front) + faoi_fn_back = (faoi_fn_back.faoi + if hasattr(faoi_fn_back, 'faoi') else faoi_fn_back) self.vf_aoi_methods = AOIMethods( faoi_fn_front, faoi_fn_back, n_integral_sections=n_aoi_integral_sections) @@ -104,19 +113,30 @@ def build_ts_vf_matrix(self, pvarray): return vf_matrix - def build_ts_vf_aoi_matrix(self, pvarray): + def build_ts_vf_aoi_matrix(self, pvarray, rho_mat): """Calculate the view factor aoi matrix elements from all PV row - surfaces to all other surfaces, only. This will not calculate + surfaces to all other surfaces, only. + If the AOI methods are available, the vf_aoi_matrix will account + for reflection losses that are AOI specific. Otherwise it will + assume that all the reflection losses are diffuse. + + Notes + ----- + When using fAOI methods, this will not calculate view factors from ground surfaces to PV row surfaces, so the users will need to run :py:meth:`~pvfactors.viewfactors.calculator.VFCalculator.build_ts_vf_matrix` first if they want the complete matrix, otherwise those entries will have zero values in them. + Parameters ---------- pvarray : :py:class:`~pvfactors.geometry.pvarray.OrderedPVArray` PV array whose timeseries view factor AOI matrix to calculate + rho_mat : np.ndarray + 2D matrix of reflectivity values for all the surfaces in the + PV array + sky. Shape = [n_ts_surfaces + 1, n_ts_surfaces + 1] Returns ------- @@ -138,16 +158,23 @@ def build_ts_vf_aoi_matrix(self, pvarray): ts_ground = pvarray.ts_ground ts_pvrows = pvarray.ts_pvrows - # Calculate vf_aoi between pvrow and ground surfaces - self.vf_aoi_methods.vf_aoi_pvrow_to_gnd(ts_pvrows, ts_ground, - tilted_to_left, - vf_aoi_matrix) - # Calculate vf_aoi between pvrows - self.vf_aoi_methods.vf_aoi_pvrow_to_pvrow(ts_pvrows, tilted_to_left, - vf_aoi_matrix) - # Calculate vf_aoi between prows and sky - self.vf_aoi_methods.vf_aoi_pvrow_to_sky(ts_pvrows, ts_ground, - tilted_to_left, vf_aoi_matrix) + if self.vf_aoi_methods is None: + # The reflection losses will be considered all diffuse. + faoi_diffuse = 1. - rho_mat + # use broadcasting + vf_aoi_matrix = faoi_diffuse * np.moveaxis(vf_aoi_matrix, -1, 0) + vf_aoi_matrix = np.moveaxis(vf_aoi_matrix, 0, -1) + else: + # Calculate vf_aoi between pvrow and ground surfaces + self.vf_aoi_methods.vf_aoi_pvrow_to_gnd(ts_pvrows, ts_ground, + tilted_to_left, + vf_aoi_matrix) + # Calculate vf_aoi between pvrows + self.vf_aoi_methods.vf_aoi_pvrow_to_pvrow( + ts_pvrows, tilted_to_left, vf_aoi_matrix) + # Calculate vf_aoi between prows and sky + self.vf_aoi_methods.vf_aoi_pvrow_to_sky( + ts_pvrows, ts_ground, tilted_to_left, vf_aoi_matrix) # Save results self.vf_aoi_matrix = vf_aoi_matrix diff --git a/pvfactors/viewfactors/vfmethods.py b/pvfactors/viewfactors/vfmethods.py index b0d17af8..59434542 100644 --- a/pvfactors/viewfactors/vfmethods.py +++ b/pvfactors/viewfactors/vfmethods.py @@ -40,9 +40,15 @@ def vf_pvrow_gnd_surf(self, ts_pvrows, ts_ground, tilted_to_left, # Front side front = ts_pvrow.front for pvrow_surf in front.all_ts_surfaces: + if pvrow_surf.is_empty: + # do no run calculation for this surface + continue ts_length = pvrow_surf.length i = pvrow_surf.index for gnd_surf in left_gnd_surfaces: + if gnd_surf.is_empty: + # do no run this calculation + continue j = gnd_surf.index vf_pvrow_to_gnd, vf_gnd_to_pvrow = ( self.vf_pvrow_surf_to_gnd_surf_obstruction_hottel( @@ -52,6 +58,9 @@ def vf_pvrow_gnd_surf(self, ts_pvrows, ts_ground, tilted_to_left, vf_matrix[i, j, :] = vf_pvrow_to_gnd vf_matrix[j, i, :] = vf_gnd_to_pvrow for gnd_surf in right_gnd_surfaces: + if gnd_surf.is_empty: + # do no run this calculation + continue j = gnd_surf.index vf_pvrow_to_gnd, vf_gnd_to_pvrow = ( self.vf_pvrow_surf_to_gnd_surf_obstruction_hottel( @@ -63,9 +72,15 @@ def vf_pvrow_gnd_surf(self, ts_pvrows, ts_ground, tilted_to_left, # Back side back = ts_pvrow.back for pvrow_surf in back.all_ts_surfaces: + if pvrow_surf.is_empty: + # do no run calculation for this surface + continue ts_length = pvrow_surf.length i = pvrow_surf.index for gnd_surf in left_gnd_surfaces: + if gnd_surf.is_empty: + # do no run this calculation + continue j = gnd_surf.index vf_pvrow_to_gnd, vf_gnd_to_pvrow = ( self.vf_pvrow_surf_to_gnd_surf_obstruction_hottel( @@ -75,6 +90,9 @@ def vf_pvrow_gnd_surf(self, ts_pvrows, ts_ground, tilted_to_left, vf_matrix[i, j, :] = vf_pvrow_to_gnd vf_matrix[j, i, :] = vf_gnd_to_pvrow for gnd_surf in right_gnd_surfaces: + if gnd_surf.is_empty: + # do no run this calculation + continue j = gnd_surf.index vf_pvrow_to_gnd, vf_gnd_to_pvrow = ( self.vf_pvrow_surf_to_gnd_surf_obstruction_hottel( @@ -185,9 +203,15 @@ def vf_pvrow_to_pvrow(self, ts_pvrows, tilted_to_left, vf_matrix): # front side front = ts_pvrow.front for surf_i in front.all_ts_surfaces: + if surf_i.is_empty: + # do no run calculation for this surface + continue i = surf_i.index length_i = surf_i.length for surf_j in right_ts_pvrow.back.all_ts_surfaces: + if surf_j.is_empty: + # do no run calculation for this surface + continue j = surf_j.index length_j = surf_j.length vf_i_to_j = self._vf_surface_to_surface( @@ -201,9 +225,15 @@ def vf_pvrow_to_pvrow(self, ts_pvrows, tilted_to_left, vf_matrix): # back side back = ts_pvrow.back for surf_i in back.all_ts_surfaces: + if surf_i.is_empty: + # do no run calculation for this surface + continue i = surf_i.index length_i = surf_i.length for surf_j in right_ts_pvrow.front.all_ts_surfaces: + if surf_j.is_empty: + # do no run calculation for this surface + continue j = surf_j.index length_j = surf_j.length vf_i_to_j = self._vf_surface_to_surface(