PV models speed up (#1212)

* updateCapacityForLifetime increase-able capacity

* update test value

* use relative tolerance for pv tests & speed up sssky_diffuse_table

* fix import of pow

* add solarpos_spa hash table

* fix test

* fix windows build

* Update save_as_JSON_test.cpp

* add lifetime table to irradproc

* update for pvsamv1

* fix pvsamv1 test results

* revert test results

* modify test results for cmod_hybrid_test

* update test results for pvwattsv8

* update PVwattsv8 tests w/ DCACRatio*

* update results for cmod_pvwattsv7 and v5

* update sam_pv_test result

* Update lib_irradproc.cpp

* minor test changes

---------

Co-authored-by: Matt Prilliman <54449384+mjprilliman@users.noreply.github.com>

shared/lib_irradproc.h

@@ -628,9 +628,68 @@ double sun_rise_and_set(double* m_rts, double* h_rts, double* delta_prime, doubl
 * \param[out] needed_values[7] zenith topocentric zenith angle (degrees)
 * \param[out] needed_values[8] azimuth topocentric azimuth angle (degrees)
 */
+
+
 void calculate_spa(double jd, double lat, double lng, double alt, double pressure, double temp,
     double delta_t, double tilt, double azm_rotation, double ascension_and_declination[2], double needed_values[9]);
 
+/**
+ * Table for storing the recently computed solarpos_spa intermediate outputs
+ * The algorithm reuses the outputs from the last 3 days or so, so the hash table is emptied every 3 days to reduce size
+ * Latitude, Longitude, Altitude, Tilt and Azimuth are not in the key because they remain constant throughout
+ */
+struct spa_table_key {
+    double jd;
+    double delta_t;
+    int pressure;
+    int temp;
+    // these are both inputs and outputs (e.g. also stored in the output vector)
+    double ascension;
+    double declination;
+
+    bool operator==(const spa_table_key &other) const
+        { return (jd == other.jd
+                && delta_t == other.delta_t
+                && pressure == other.pressure
+                && temp == other.temp
+                && ascension == other.ascension
+                && declination == other.declination
+                );
+        }
+
+    spa_table_key(double j, double dt, double p, double t, double a, double d):
+        jd(j), delta_t(dt), ascension(a), declination(d)
+    {
+        int pressure_bucket = 10;
+        pressure = ((int)(p + pressure_bucket/2) / pressure_bucket) * pressure_bucket;
+        pressure = (int)pressure;
+        int temp_bucket = 5;
+        temp = ((int)(t + temp_bucket / 2) / temp_bucket) * temp_bucket;
+        temp = (int)temp;
+    }
+};
+
+template <>
+struct std::hash<spa_table_key>
+{
+  std::size_t operator()(const spa_table_key& k) const
+  {
+    using std::hash;
+    // Compute individual hash values for first, second, etc and combine them using XOR and bit shifting:
+    return 
+    ((((
+            ((((hash<double>()(k.jd)
+             ^ (hash<double>()(k.delta_t) << 1)) >> 1)
+             ^ (hash<int>()(k.pressure) << 1)) >> 1)
+             ^ (hash<int>()(k.temp) << 1)) >> 1)
+             ^ (hash<double>()(k.ascension) << 1)) >> 1)
+             ^ (hash<double>()(k.declination) << 1)
+             ;
+  }
+};
+
+void clear_spa_table();
+
 /**
 *
 *	calculate_eot_and_sun_rise_transit_set function calculates Equation of Time (EoT) and sunrise and
@@ -1048,6 +1107,8 @@ class irrad
     std::vector<double> planeOfArrayIrradianceRearSpatialCS;  ///< Spatial rear side clearsky plane-of-array irradiance (W/m2), where index 0 is at row bottom
     std::vector<double> groundIrradianceSpatial;            ///< Spatial irradiance incident on the ground in between rows, where index 0 is towards front of array
 
+    std::vector<std::vector<double>> solarpos_outputs_for_lifetime;     ///< Table of solarpos outputs stored for lifetime simulations
+    void storeSolarposOutputs();
 public:
 
     /// Directive to indicate that if delt_hr is less than zero, do not interpolate sunrise and sunset hours
@@ -1063,13 +1124,11 @@ class irrad
     static const int groundIrradOutputRes = 10;
 
     /// Default class constructor, calls setup()
-    irrad(weather_record wr, weather_header wh,
+    irrad(weather_header wh,
         int skyModel, int radiationModeIn, int trackModeIn,
-        bool useWeatherFileAlbedo, bool instantaneousWeather, bool backtrackingEnabled, bool forceToStowIn,
+        bool instantaneousWeather, bool backtrackingEnabled, bool forceToStowIn,
         double dtHour, double tiltDegrees, double azimuthDegrees, double trackerRotationLimitDegrees, double stowAngleDegreesIn,
-        double groundCoverageRatio, double slopeTilt, double slopeAzm, std::vector<double> monthlyTiltDegrees, std::vector<double> userSpecifiedAlbedo,
-        poaDecompReq* poaAllIn,
-        bool useSpatialAlbedos = false, const util::matrix_t<double>* userSpecifiedSpatialAlbedos = nullptr, bool enableSubhourlyClipping = false, bool useCustomRotAngles = false, double customRotAngle = 0);
+        double groundCoverageRatio, double slopeTilt, double slopeAzm, poaDecompReq *poaAllIn, bool enableSubhourlyClipping = false);
 
     /// Construct the irrad class with an Irradiance_IO() object and Subarray_IO() object
     irrad();
@@ -1080,6 +1139,9 @@ class irrad
     /// Validation method to verify member data is within acceptable ranges
     int check();
 
+    /// Initialize solarpos_outputs_lifetime for lifetime simulations
+    void setup_solarpos_outputs_for_lifetime(size_t n_steps_per_year);
+
     /// Set the time for the irradiance processor
     void set_time(int year, int month, int day, int hour, double minute, double delt_hr);
 
@@ -1118,6 +1180,11 @@ class irrad
     /// Function to overwrite internally calculated sun position values, primarily to enable testing against other libraries using different sun position calculations
     void set_sun_component(size_t index, double value);
 
+    /// Function to set time, albedo, tracking and optional inputs from weather record
+    void set_from_weather_record(weather_record wf, weather_header hdr, int trackModeIn, std::vector<double>& monthlyTiltDegrees, 
+            bool useWeatherFileAlbedo, std::vector<double>& userSpecifiedAlbedo, poaDecompReq *poaAllIn, bool useSpatialAlbedos, const util::matrix_t<double>* userSpecifiedSpatialAlbedos, 
+            bool useCustomRotAngles = false, double customRotAngle = 0);
+
     /// Run the irradiance processor and calculate the plane-of-array irradiance and diffuse components of irradiance
     int calc();
 
@@ -1216,6 +1283,9 @@ class irrad
     /// Return the front surface irradiances, used by \link calc_rear_side()
     void getFrontSurfaceIrradiances(double pvBackShadeFraction, double rowToRow, double verticalHeight, double clearanceGround, double distanceBetweenRows, double horizontalLength, std::vector<double> frontGroundGHI, std::vector<double>& frontIrradiance, double& frontAverageIrradiance, std::vector<double>& frontReflected);
 
+    /// Return the solarpos outputs for a given timestep
+    bool getStoredSolarposOutputs();
+
     enum RADMODE { DN_DF, DN_GH, GH_DF, POA_R, POA_P };
     enum SKYMODEL { ISOTROPIC, HDKR, PEREZ };
     enum TRACKING { FIXED_TILT, SINGLE_AXIS, TWO_AXIS, AZIMUTH_AXIS, SEASONAL_TILT };

shared/lib_pvshade.cpp

@@ -291,27 +291,28 @@ void selfshade_xs_horstr(bool landscape,
 
 // Accessor for sky diffuse derates for the given surface_tilt. If the value doesn't exist in the table, it is computed.
 double sssky_diffuse_table::lookup(double surface_tilt) {
-    char buf[8];
-    sprintf(buf, "%.3f", surface_tilt);
-    if (derates_table.find(buf) != derates_table.end())
-        return derates_table[buf];
+    int surface_tilt_key = int(surface_tilt * derates_table_digits_multiplier);
+
+    if (derates_table.find(surface_tilt_key) != derates_table.end())
+        return derates_table[surface_tilt_key];
     return compute(surface_tilt);
 }
 
 double sssky_diffuse_table::compute(double surface_tilt) {
     if (gcr == 0)
         throw std::runtime_error("sssky_diffuse_table::compute error: gcr required in initialization");
     // sky diffuse reduction
-    double step = 1.0 / 1000.0;
+	const size_t n_steps = 250;
+    double step = 1.0 / (double)n_steps;
     double skydiff = 0.0;
     double tand_stilt = tand(surface_tilt);
     double sind_stilt = sind(surface_tilt);
     double Asky = M_PI + M_PI / pow((1 + pow(tand_stilt, 2)), 0.5);
-    double arg[1000];
-    double gamma[1000];
-    double tan_tilt_gamma[1000];
-    double Asky_shade[1000];
-    for (int n = 0; n < 1000; n++)
+    double arg[n_steps];
+    double gamma[n_steps];
+    double tan_tilt_gamma[n_steps];
+    double Asky_shade[n_steps];
+    for (int n = 0; n < n_steps; n++)
     {
         if (surface_tilt != 0)
             arg[n] = (1 / tand_stilt) - (1 / (gcr * sind_stilt * (1 - n * step)));
@@ -331,9 +332,8 @@ double sssky_diffuse_table::compute(double surface_tilt) {
         else {}
         skydiff += (Asky_shade[n] / Asky) * step;
     }
-    char buf[8];
-    sprintf(buf, "%.3f", surface_tilt);
-    derates_table[buf] = skydiff;
+    int surface_tilt_key = int(surface_tilt * derates_table_digits_multiplier);
+	derates_table[surface_tilt_key] = skydiff;
     return skydiff;
 }
 

shared/lib_pvshade.h

@@ -36,6 +36,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include <string>
 #include <unordered_map>
+#include <cmath>
 #include "lib_util.h"
 
 
@@ -73,7 +74,9 @@ struct ssoutputs	// self-shading outputs
 // added to removing duplicate computations for speed up (https://github.com/NREL/ssc/issues/384)
 class sssky_diffuse_table
 {
-    std::unordered_map<std::string, double> derates_table;      // stores pairs of surface tilts and derates
+    std::unordered_map<int, double> derates_table;      // stores pairs of surface tilts and derates
+	size_t derates_table_digits = 0;
+	size_t derates_table_digits_multiplier = std::pow(10, derates_table_digits);
     double gcr;                                                 // 0.01 - 0.99
 
     double compute(double surface_tilt);

shared/lib_util.cpp

@@ -1255,6 +1255,13 @@ size_t util::yearOneIndex(double dtHour, size_t lifetimeIndex)
 	return indexYearOne;
 }
 
+size_t util::yearIndex(size_t year, size_t month, size_t day, size_t hour, double minute, size_t steps_per_hour)
+{
+	size_t mins_in_step = (size_t)(60 / steps_per_hour);
+    size_t step_of_hour = (size_t)(minute / mins_in_step);
+    return util::lifetimeIndex((size_t)year, (size_t)util::hour_of_year(month, day, hour), step_of_hour, steps_per_hour);
+}
+
 std::vector<double> util::frequency_table(double* values, size_t n_vals, double bin_width)
 {
     if (!values)

shared/lib_util.h

@@ -143,6 +143,7 @@ namespace util
 	bool weekday(size_t hour_of_year); /* return true if is a weekday, assuming first hour of year is Monday at 12 am*/
 	size_t lifetimeIndex(size_t year, size_t hour_of_year, size_t step_of_hour, size_t steps_per_hour);
 	size_t yearOneIndex(double dtHour, size_t lifetimeIndex);
+	size_t yearIndex(size_t year, size_t month, size_t day, size_t hour, double minute, size_t step_per_hour);
 
 	int schedule_char_to_int( char c );
 	std::string schedule_int_to_month( int m );

shared/lib_weatherfile.cpp

@@ -1542,6 +1542,8 @@ bool weatherfile::open(const std::string& file, bool header_only)
         }
     }
 
+    m_startYear = m_columns[YEAR].data[0];
+
     if (!header_only) {
         start_hours_at_0();
 

shared/lib_weatherfile.h

@@ -201,6 +201,8 @@ class weather_data_provider
 	double tz() { return header().tz; }
 	double elev() { return header().elev; }
 
+	double start_year() {return m_startYear; }
+
     bool check_hour_of_year(int hour, int line);
 
 	// virtual functions specific to weather data source