Merge pull request #120 from respec/develop

Enhancements for custom python SPECL, small GENER fixes

.gitignore

@@ -62,3 +62,6 @@ tests/GLWACSO/HSP2results/hspp007.hdf
 tests/GLWACSO/HSPFresults/hspf006.HBN
 tests/GLWACSO/HSP2results/hspp007.uci
 tests/test_report_conversion.html
+tests/land_spec/hwmA51800.h5
+tests/testcbp/HSP2results/PL3_5250_0001.h5
+tests/testcbp/HSP2results/*.csv

HSP2/ADCALC.py

@@ -110,6 +110,9 @@ def _adcalc_(ui, ts):
 			ts['EOVOL' + str(index + 1)] = EOVOL[:, index]
 
 	ROS = ui['ROS']
+	OS  = zeros(nexits)
+	for index in range(nexits):
+		OS[index] = ui['OS' + str(index + 1)]
 
 	# external time series
 	O = zeros((simlen, nexits))
@@ -123,7 +126,7 @@ def _adcalc_(ui, ts):
 		vols = VOL[loop-1] * VFACT  if loop > 0 else vol
 
 		o  = O[loop]
-		os = O[loop-1] if loop > 0 else array([ROS])
+		os = O[loop-1] if loop > 0 else OS
 		ro = 0.0
 		ros= 0.0
 		for index in range(nexits):

HSP2/GENER.py

@@ -1,14 +1,15 @@
 import numpy as np
 import pandas as pd
-from typing import Dict
+from numpy import zeros
+from numba.typed import Dict
 
 class Gener():
     """
     Partial implementation of the GENER module. 
-    Currently supports OPCODES 1-7, 9-23, & 25-26
+    Currently supports OPCODES 1-7, 9-26
     """
 
-    def __init__(self, segment: str, copies: Dict, geners: Dict, ddlinks: Dict, ddgener: Dict) -> None:
+    def __init__(self, segment: str, siminfo: Dict, copies: Dict, geners: Dict, ddlinks: Dict, ddgener: Dict) -> None:
         self.ts_input_1 = pd.Series() # type: pd.Series
         self. ts_input_2 = pd.Series() # type: pd.Series
         self.ts_output = pd.Series()  # type: pd.Series
@@ -19,16 +20,26 @@ def __init__(self, segment: str, copies: Dict, geners: Dict, ddlinks: Dict, ddge
         if self.opcode in [9,10,11,24,25,26]:
             self.k = ddgener['PARM'][segment]
 
+        # special case for k as constant case 24
+        if self.opcode == 24:
+            start, stop, steps = siminfo['start'], siminfo['stop'], siminfo['steps']
+            ts = zeros(steps)
+            self.ts_input_1 = ts
+
         for link in ddlinks[segment]:
             ts = pd.Series()
             if link.SVOL == 'COPY': 
                 copy = copies[link.SVOLNO]
                 ts = copy.get_ts(link.SMEMN,link.SMEMSB1)
                 if link.MFACTOR != 1: ts *= link.MFACTOR
             elif link.SVOL == 'GENER':
-                gener = geners[link.SVOLNO]
-                ts = gener.get_ts()
-                if link.MFACTOR != 1: ts *= link.MFACTOR
+                if link.SVOLNO in geners:
+                    gener = geners[link.SVOLNO]
+                    ts = gener.get_ts()
+                    if link.MFACTOR != 1: ts *= link.MFACTOR
+                else:
+                    raise NotImplementedError(
+                        f"Invalid SVOL. This GENER operation does not exist. '{link.SVOLNO}'")
             else:
                 raise NotImplementedError(f"Invalid SVOL. GENER module does not currently support reading TimeSeries for '{link.SVOL}'")
 
@@ -135,9 +146,8 @@ def _opcode23(self) -> pd.Series:
         return ts_out
 
     def _opcode24(self) -> pd.Series:
-        #skip for now 
-        #would need to figure out timeseries length component
-        raise NotImplementedError("GENER OPCODE 24 is not currently supported")
+        # ts_input_1 is all zeros, just add k to it
+        return np.add(self.ts_input_1, self.k)
 
     def _opcode25(self) -> pd.Series:
         return np.maximum(self.ts_input_1, self.k)

HSP2/HTRCH.py

@@ -287,7 +287,7 @@ def _htrch_(ui, ts):
 		# calculate heat transfer rates for water surface; units are kcal/m2.ivl
 
 		# get quantity of precipitation and convert ft/ivl to m/ivl,
-		prec = PREC[loop] / 12.0
+		prec = PREC[loop] * 0.0833  # / 12.0 to match HSPF precision
 		if prec > 0.0:
 			if uunits == 1:
 				mprec = prec /3.2808

HSP2/HYDR.py

@@ -12,12 +12,14 @@
 '''
 
 
-from numpy import zeros, any, full, nan, array, int64
+from numpy import zeros, any, full, nan, array, int64, arange
 from pandas import DataFrame
 from math import sqrt, log10
 from numba import njit
 from numba.typed import List
 from HSP2.utilities import initm, make_numba_dict
+from HSP2.state import *
+from HSP2.SPECL import specl
 
 
 ERRMSGS =('HYDR: SOLVE equations are indeterminate',             #ERRMSG0
@@ -30,15 +32,18 @@
 MAXLOOPS  = 100     # newton method exit tolerance
 
 
-def hydr(io_manager, siminfo, uci, ts, ftables):
+def hydr(io_manager, siminfo, uci, ts, ftables, state):
     ''' find the state of the reach/reservoir at the end of the time interval
     and the outflows during the interval
 
-    CALL: hydr(store, general, ui, ts)
+    CALL: hydr(store, general, ui, ts, specactions)
        store is the Pandas/PyTable open store
        general is a dictionary with simulation level infor (OP_SEQUENCE for example)
        ui is a dictionary with RID specific HSPF UCI like data
-       ts is a dictionary with RID specific timeseries'''
+       ts is a dictionary with RID specific timeseries
+       state is a dictionary that contains all dynamic code dictionaries such as: 
+       - specactions is a dictionary with all special actions
+    '''
 
     steps   = siminfo['steps']                # number of simulation points
     uunits  = siminfo['units']
@@ -117,22 +122,48 @@ def hydr(io_manager, siminfo, uci, ts, ftables):
     for i in range(nexits):
         Olabels.append(f'O{i+1}')
         OVOLlabels.append(f'OVOL{i+1}')
-
+
+    # state_info is some generic things about the simulation 
+    state_info = Dict.empty(key_type=types.unicode_type, value_type=types.unicode_type)
+    state_info['operation'], state_info['segment'], state_info['activity'] = state['operation'], state['segment'], state['activity']
+    state_info['domain'], state_info['state_step_hydr'] = state['domain'], state['state_step_hydr']
+    hsp2_local_py = state['hsp2_local_py']
+    # It appears necessary to load this here, instead of from main.py, otherwise,
+    # _hydr_() does not recognize the function state_step_hydr()? 
+    if (hsp2_local_py != False):
+        from hsp2_local_py import state_step_hydr
+    else:
+        from HSP2.state_fn_defaults import state_step_hydr
+    # must split dicts out of state Dict since numba cannot handle mixed-type nested Dicts
+    state_ix, dict_ix, ts_ix = state['state_ix'], state['dict_ix'], state['ts_ix']
+    state_paths = state['state_paths']
+    # initialize the hydr paths in case they don't already reside here
+    hydr_init_ix(state_ix, state_paths, state['domain'])
+
+    ###########################################################################
+    # specactions - special actions code TBD
+    ###########################################################################
+    specactions = make_numba_dict(state['specactions']) # Note: all values coverted to float automatically
+
+    ###########################################################################
+    # Do the simulation with _hydr_()
     ###########################################################################
-    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels)                  # run reaches simulation code
+    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, specactions, state_step_hydr) # run reaches simulation code
     ###########################################################################
 
     if 'O'    in ts:  del ts['O']
     if 'OVOL' in ts:  del ts['OVOL']
 
     # save initial outflow(s) from reach:
     uci['PARAMETERS']['ROS'] = ui['ROS']
+    for i in range(nexits):
+        uci['PARAMETERS']['OS'+str(i+1)] = ui['OS'+str(i+1)]
 
     return errors, ERRMSGS
 
 
 @njit(cache=True)
-def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels):
+def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, specactions, state_step_hydr):
     errors = zeros(int(ui['errlen'])).astype(int64)
 
     steps  = int(ui['steps'])            # number of simulation steps
@@ -255,15 +286,54 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels):
 
     # store initial outflow from reach:
     ui['ROS'] = ro
-
+    for index in range(nexits):
+        ui['OS' + str(index + 1)] = o[index]
+
+    # other initial vars
+    rovol = 0.0
+    volev = 0.0
+    IVOL0   = ts['IVOL']                   # the actual inflow in simulation native units 
+    # prepare for dynamic state
+    hydr_ix = hydr_get_ix(state_ix, state_paths, state_info['domain'])
+    # these are integer placeholders faster than calling the array look each timestep
+    o1_ix, o2_ix, o3_ix, ivol_ix = hydr_ix['O1'], hydr_ix['O2'], hydr_ix['O3'], hydr_ix['IVOL']
+    ro_ix, rovol_ix, volev_ix, vol_ix = hydr_ix['RO'], hydr_ix['ROVOL'], hydr_ix['VOLEV'], hydr_ix['VOL']
+    # handle varying length outdgt
+    out_ix = arange(nexits)
+    if nexits > 0:
+        out_ix[0] = o1_ix
+    if nexits > 1:
+        out_ix[1] = o2_ix
+    if nexits > 2:
+        out_ix[2] = o3_ix
     # HYDR (except where noted)
     for step in range(steps):
+        # call specl
+        specl(ui, ts, step, state_info, state_paths, state_ix, specactions)
         convf  = CONVF[step]
         outdgt[:] = OUTDGT[step, :]
         colind[:] = COLIND[step, :]
         roseff = ro
         oseff[:] = o[:]
-
+        # set state_ix with value of local state variables and/or needed vars
+        # Note: we pass IVOL0, not IVOL here since IVOL has been converted to different units
+        state_ix[ro_ix], state_ix[rovol_ix] = ro, rovol
+        di = 0
+        for oi in range(nexits):
+            state_ix[out_ix[oi]] = outdgt[oi] 
+        state_ix[vol_ix], state_ix[ivol_ix] = vol, IVOL0[step]
+        state_ix[volev_ix] = volev
+        # Execute dynamic code if enabled
+        if (state_info['state_step_hydr'] == 'enabled'):
+            state_step_hydr(state_info, state_paths, state_ix, dict_ix, ts_ix, hydr_ix, step)
+            # Do write-backs for editable STATE variables
+            # OUTDGT is writeable
+            for oi in range(nexits):
+                outdgt[oi] = state_ix[out_ix[oi]]
+            # IVOL is writeable. 
+            # Note: we must convert IVOL to the units expected in _hydr_ 
+            # maybe routines should do this, and this is not needed (but pass VFACT in state)
+            IVOL[step] = state_ix[ivol_ix] * VFACT
         # vols, sas variables and their initializations  not needed.
         if irexit >= 0:             # irrigation exit is set, zero based number
             if rirwdl > 0.0:  # equivalent to OVOL for the irrigation exit
@@ -602,4 +672,5 @@ def expand_HYDR_masslinks(flags, uci, dat, recs):
         rec['TMEMSB2'] = dat.TMEMSB2
         rec['SVOL'] = dat.SVOL
         recs.append(rec)
-    return recs
+    return recs
+

HSP2/IQUAL.py

@@ -7,7 +7,7 @@
 from math import exp
 from numpy import zeros, where, full, float64, int64
 from numba import njit
-from HSP2.utilities import initm, make_numba_dict, hourflag
+from HSP2.utilities import initm, make_numba_dict, hourflag, initmdiv
 
 
 ''' DESIGN NOTES
@@ -30,7 +30,7 @@ def iqual(io_manager, siminfo, uci, ts):
 	nquals = 1
 	if 'PARAMETERS' in uci:
 		if 'NQUAL' in uci['PARAMETERS']:
-			nquals = uci['PARAMETERS']['NQUAL']
+			nquals = int(uci['PARAMETERS']['NQUAL'])
 	constituents = []
 	for index in range(nquals):
 		iqual = str(index + 1)
@@ -70,6 +70,10 @@ def iqual(io_manager, siminfo, uci, ts):
 		ts['ACQOP' + str(index)] = initm(siminfo, uci, ui_flags['VQOFG'], 'IQUAL' + str(index) + '_MONTHLY/ACQOP', ui_parms['ACQOP'])
 		ts['SQOLIM' + str(index)] = initm(siminfo, uci, ui_flags['VQOFG'], 'IQUAL' + str(index) + '_MONTHLY/SQOLIM', ui_parms['SQOLIM'])
 
+		ts['REMQOP' + str(index)] = initmdiv(siminfo, uci, ui_flags['VQOFG'], 'IQUAL' + str(index) + '_MONTHLY/ACQOP',
+											 'IQUAL' + str(index) + '_MONTHLY/SQOLIM', ui_parms['ACQOP'],
+											 ui_parms['SQOLIM'])
+
 		iqadfgf = 0
 		iqadfgc = 0
 		ts['IQADFX' + str(index)] = zeros(simlen)
@@ -169,6 +173,7 @@ def _iqual_(ui, ts):
 		POTFW  = ts['POTFW'  + str(index)]
 		ACQOP  = ts['ACQOP'  + str(index)]
 		SQOLIM = ts['SQOLIM' + str(index)]
+		REMQOP = ts['REMQOP' + str(index)]
 		IQADFX = ts['IQADFX' + str(index)]
 		IQADCN = ts['IQADCN' + str(index)]
 
@@ -218,7 +223,7 @@ def _iqual_(ui, ts):
 					# washof ()
 					''' Simulate accumulation of a quality constituent on the land surface and its removal using a constant unit rate and by direct washoff by overland flow'''
 					if dayfg == 1:
-						remqop = acqop / SQOLIM[loop]
+						remqop = REMQOP[loop]
 						if QSOFG == 1 :   #update storage due to accumulation and removal which occurs independent of runoff - units are qty/acre
 							sqo = acqop + sqo * (1.0 - remqop)
 

HSP2/PQUAL.py

@@ -7,7 +7,7 @@
 from math import exp
 from numpy import zeros, where, full, float64, int64
 from numba import njit
-from HSP2.utilities import initm, make_numba_dict, hourflag
+from HSP2.utilities import initm, make_numba_dict, hourflag, initmdiv
 
 ''' DESIGN NOTES
 Each constituent will be in its own subdirectory in the HDF5 file.
@@ -33,7 +33,7 @@ def pqual(io_manager, siminfo, uci, ts):
 	nquals = 1
 	if 'PARAMETERS' in uci:
 		if 'NQUAL' in uci['PARAMETERS']:
-			nquals = uci['PARAMETERS']['NQUAL']
+			nquals = int(uci['PARAMETERS']['NQUAL'])
 	constituents = []
 	for index in range(nquals):
 		pqual = str(index + 1)
@@ -81,6 +81,10 @@ def pqual(io_manager, siminfo, uci, ts):
 		ts['IOQCP' + str(index)] = initm(siminfo, uci, ui_flags['VIQCFG'], 'PQUAL' + str(index) + '_MONTHLY/IOQC', ui_parms['IOQC'])
 		ts['AOQCP' + str(index)] = initm(siminfo, uci, ui_flags['VAQCFG'], 'PQUAL' + str(index) + '_MONTHLY/AOQC', ui_parms['AOQC'])
 
+		ts['REMQOP' + str(index)] = initmdiv(siminfo, uci, ui_flags['VQOFG'], 'PQUAL' + str(index) + '_MONTHLY/ACQOP',
+											 'PQUAL' + str(index) + '_MONTHLY/SQOLIM', ui_parms['ACQOP'],
+											 ui_parms['SQOLIM'])
+
 		pqadfgf = 0
 		pqadfgc = 0
 		ts['PQADFX' + str(index)] = zeros(simlen)
@@ -204,6 +208,7 @@ def _pqual_(ui, ts):
 		POTFS  = ts['POTFS' + str(index)]
 		ACQOP  = ts['ACQOP' + str(index)]
 		SQOLIM = ts['SQOLIM' + str(index)]
+		REMQOP = ts['REMQOP' + str(index)]
 		IOQCP  = ts['IOQCP' + str(index)]
 		AOQCP  = ts['AOQCP' + str(index)]
 		PQADFX = ts['PQADFX' + str(index)]
@@ -226,6 +231,7 @@ def _pqual_(ui, ts):
 			potfs  = POTFS[loop]
 			acqop  = ACQOP[loop]
 			sqolim = SQOLIM[loop]
+			remqop = REMQOP[loop]
 			ioqc   = IOQCP[loop]
 			aoqc   = AOQCP[loop]
 			iliqc  = ILIQC[loop]
@@ -272,7 +278,7 @@ def _pqual_(ui, ts):
 				# qualof()
 				''' Simulate accumulation of a quality constituent on the land surface and its removal by a constant unit rate and by overland flow'''	
 				if dayfg:
-					remqop = acqop / sqolim
+					# remqop = acqop / sqolim
 
 					if QSOFG == 1:
 						# update storage due to accumulation and removal which occurs independent of runoff - units are qty/acre

HSP2/PSTEMP.py

@@ -119,8 +119,8 @@ def _pstemp_(ui, ts):
 
 	if uunits != 2:
 		if TSOPFG == 1:
-			ULTP1= (ULTP1 - 32.0) * 0.555
-			LGTP1= (LGTP1 - 32.0) * 0.555
+			ULTP1= (ULTP1 - 32.0) * 0.5555 # trying to match HSPF precision here
+			LGTP1= (LGTP1 - 32.0) * 0.5555
 		else:
 			ULTP2=  0.555 * ULTP2
 			LGTP2 = 0.555 * LGTP2

HSP2/PWATER.py

@@ -516,7 +516,7 @@ def _pwater_(ui, ts):
         elif agws > 1.0e-20:
             agwo = kgw * agws  # enough water to have outflow
 
-        if agwo < 0.0:
+        if agwo < 1.0e-12:   # beyond simgle precision math, HSPF will have zero
             agwo = 0.0
 
         # no remaining water - this should happen only with hwtfg=1 it may
@@ -595,9 +595,9 @@ def _pwater_(ui, ts):
 
             if abs(kvary) > 0.0:
                 gwvs -= agwet   # update variable storage
-                if gwvs < -0.02:
-                    errors[5] += 1.0   # ERRMSG5: GWVS < -0.02, set to zero
-                    gwvs = 0.0
+                # if gwvs < -0.02:   # this check is commented out in HSPF v12+
+                #     errors[5] += 1.0   # ERRMSG5: GWVS < -0.02, set to zero
+                #     gwvs = 0.0
             taet   += agwet
             rempet -= agwet