junction.py

Go to the documentation of this file.

# -*- coding: utf-8 -*-
"""
Created on Wed Apr 29 22:04:46 2020
 
@author: cd
"""
 
 
import abc
import numpy as np
from efo.general import Constants
from efo.time import TimeBase
from efo.hypso import HypsoBase
from efo.outlet import OutletRatedElev
from efo.rule_compliance import RuleComplianceBase
 
 
class JunctionBase(metaclass=abc.ABCMeta):
    def __init__(self, name, time, qIn=[]):
        self.namename = name
        super().__init__()
        # Set class properties
        self.qInqIn = []
        for qInCur in qIn:
            self.append_qinappend_qin(qInCur)
        # TODO: Should use error handling instead
        if issubclass(type(time), TimeBase):
            self.TT = time
            self.qOutqOut = np.full(self.TT.nSteps, np.nan)
            self.qOutqOut[0] = 0.
            self.qInTotqInTot = np.full(self.TT.nSteps, np.nan)
    
    def append_qin(self, qIn):
        # TODO: Use error handling to check Qin class
        self.qInqIn.append(qIn)
 
    def _get_qin(self, *, tsOffset=0):
        qTot = 0.
        for qInCur in self.qInqIn:
            qTot += qInCur.get_qin(tsOffset)
        self.qInTotqInTot[self.TT.step] = qTot
        return qTot
    
    def get_qout(self, *, tsOffset=0):
        return self.qOutqOut[min(self.TT.end, max(0, self.TT.step + tsOffset))].item()
 
    @abc.abstractmethod
    def calc_qout(self):
        qOut = max(0., self._get_qin_get_qin())
        # TODO: Can qOut be a property to control writing access?
        self.qOutqOut[self.TT.step] = qOut
        return qOut
    
    def set_qout(self, qOutSpecified, *, tsOffset=0):
        self.qOutqOut[min(self.TT.end, self.TT.step + tsOffset)] = qOutSpecified
 
    @classmethod
    def __subclasshook__(cls, C):
        if cls is JunctionBase:
            attrs = set(dir(C))
            if set(cls.__abstractmethods__) <= attrs:
                return True
        return NotImplemented
 
 
class Junction(JunctionBase):
    def __init__(self, name, time, qIn=[]):
        # Call super class constructor
        super().__init__(name, time, qIn=qIn)
        
    def calc_qout(self):
        return super().calc_qout()
        
        
class JunctionRegulated(JunctionBase):
    def __init__(self, name, time, qIn,
                 ruleMinQ=None, ruleMaxQ=None, ruleDiversion=[]):
        # Call super class constructor
        super().__init__(name, time, qIn)
        # Set class properties
        # Need to make sure these are max and min rules and diversion rules
        # TODO: Should this also have a qMax and qMin like ReservoirJunction?
        # TODO: Instead of having a ruleMinQ and ruleMaxQ could this also have...
        # ... a rule stack
        self.ruleMinQruleMinQ = ruleMinQ
        self.ruleMaxQruleMaxQ = ruleMaxQ
        self.ruleDiversionruleDiversion = ruleDiversion
        if ruleDiversion: self.qDivqDiv = np.full(self.TT.nSteps, np.nan)
        self.continuitycontinuity = np.full(self.TT.nSteps, np.nan)
 
    # TODO: Property?
    def set_rule_minq(self, ruleMinQ):
        self.ruleMinQruleMinQ = ruleMinQ
 
    # TODO: Property?
    def set_rule_maxq(self, ruleMaxQ):
        self.ruleMaxQruleMaxQ = ruleMaxQ
    
    def _calc_q(self, qIn, qDiv, rule):
        qOut = qIn - qDiv
        if qOut < 0.:
            qDiv += qOut
            rule.release[self.TT.step] = qDiv
            qOut = 0.
        return qOut, qDiv
   
    def calc_delta(self, *, ruleType=RuleComplianceBase.MAX, tsOffset=0):
        if ruleType == RuleComplianceBase.MAX and self.ruleMaxQruleMaxQ:
            targetQ = self.ruleMaxQruleMaxQ.get_qcomp(tsOffset=tsOffset)
        elif ruleType == RuleComplianceBase.MIN and self.ruleMinQruleMinQ:
            targetQ = self.ruleMinQruleMinQ.get_qcomp(tsOffset=tsOffset)
        else:
            targetQ = np.nan
        qDiv = 0.
        qInInit = self._get_qin_get_qin(tsOffset=tsOffset)
        if self.ruleDiversionruleDiversion:
            for curRule in self.ruleDiversionruleDiversion:
                qDiv += curRule.get_qcomp(qIn=qInInit, tsOffset=tsOffset)
        qInNet = qInInit - qDiv
        self.qOutqOut[min(self.TT.end, self.TT.step + tsOffset)] = max(0., qInNet)
        return qInNet - targetQ if ruleType == RuleComplianceBase.MAX else targetQ - qInNet
 
    def calc_qout(self):
        qInTot = super().calc_qout()
        # self.qInNet[self.T.step] = qInTot
        qDiv = 0.
        qOut = qInTot
        qMin = self.ruleMinQruleMinQ.get_qcomp(qIn=qInTot) if self.ruleMinQruleMinQ is not None else 0.
        qDivMax = qInTot - qMin
        if self.ruleDiversionruleDiversion:
            for curRule in self.ruleDiversionruleDiversion:
                curDiv = curRule.calc_release(rlsProposed=max(0., qDivMax-qDiv), qIn=qInTot)
                qOut, curDiv = self._calc_q_calc_q(qOut, curDiv, curRule)
                qDiv += curDiv
        # if self.ruleDiversion: self.qDiv[self.T.step] = qDiv
            self.qDivqDiv[self.TT.step] = qDiv
        # qOut = qInTot - qDiv
        self.continuitycontinuity[self.TT.step] = qOut + qDiv - qInTot
        self.qOutqOut[self.TT.step] = qOut
        return qOut
 
 
 
# TODO: Maybe make the ruleStack it's own class that can have a method to evaluate the rules...
# ...this could also be used to map the rules by priority.
class ReservoirJunction(JunctionRegulated):
    def __init__(self, name, time, qIn, storInit, constants, 
                 *,hypso=None, ruleStack=[], evapObj=None, seepageObj=None, 
                 outletCtrl=None, outletUnctrl=None, ruleEmgc=None, storMax=None):
        # Call super class constructor
        super().__init__(name, time, qIn)
        # Set class properties
        self.ruleStackruleStack = ruleStack
        self.set_ctrl_outletset_ctrl_outlet(outletCtrl)
        self.set_unctrl_outletset_unctrl_outlet(outletUnctrl)
        self.ruleEmgcruleEmgc = ruleEmgc
        # TODO: This could be a guide curve rule or just a scalar
        self.storMaxstorMax = storMax
        # TODO: Should use error handling instead
        if issubclass(type(constants), Constants):
            self.constantsconstants = constants
        if issubclass(type(hypso), HypsoBase):
            self.hypsohypso = hypso
        # Create storage array
        self.storstor = np.full(self.TT.nSteps, np.nan)
        self.storstor[0] = storInit
        # Other propoerties
        self.set_evapset_evap(evapObj)
        self.set_seepageset_seepage(seepageObj)
        # self.continuity = np.full(self.T.nSteps, np.nan)
        # self.qInNet = np.full(self.T.nSteps, np.nan)
 
    def create_ctrl_outlet(self, name, elev, qOutlet):
        # TODO: Use error handling to check values passed
        outletCtrl = OutletRatedElev(name, self.TT, self.hypsohypso, elev, qOutlet)
        self.set_ctrl_outletset_ctrl_outlet(outletCtrl)
 
    def create_unctrl_outlet(self, name, elev, qOutlet):
        # TODO: Use error handling to check values passed
        outletUnCtrl = OutletRatedElev(name, self.TT, self.hypsohypso, elev, qOutlet)
        self.set_unctrl_outletset_unctrl_outlet(outletUnCtrl)
        
    def set_ctrl_outlet(self, outletCtrl):
        # TODO: Use error handling to check values passed
        self.outletCtrloutletCtrl = outletCtrl
        if outletCtrl: 
            self.rlsCtrlrlsCtrl = np.full(self.TT.nSteps, np.nan)
            self.rlsCtrlrlsCtrl[0] = 0.
            self.rlsMinrlsMin = np.full(self.TT.nSteps, np.nan)
            self.rlsMaxrlsMax = np.full(self.TT.nSteps, np.nan)
 
    def set_unctrl_outlet(self, outletUnctrl):
        # TODO: Use error handling to check values passed
        self.outletUnCtrloutletUnCtrl = outletUnctrl
        if outletUnctrl: self.rlsUnCtrlrlsUnCtrl = np.full(self.TT.nSteps, np.nan)
 
    def set_evap(self, evapObj):
        # TODO: Use error handling to check values passed
        self.evapObjevapObj = evapObj
        if evapObj: self.lossEvaplossEvap = np.full(self.TT.nSteps, np.nan)
        
    def set_seepage(self, seepageObj):
        # TODO: Use error handling to check values passed
        self.seepageObjseepageObj = seepageObj
        if seepageObj: self.lossSeeplossSeep = np.full(self.TT.nSteps, np.nan)
 
    def set_rule_emgc(self, ruleEmgc):
        # TODO: Add error handling
        self.ruleEmgcruleEmgc = ruleEmgc
        
    def set_stor_init(self, storInit):
        self.storstor[0] = storInit
        
    def set_ruleStack(self, ruleStack):
        # TODO: Error handling needed here
        if type(ruleStack) == list:
            self.ruleStackruleStack = ruleStack
 
    def append_rule(self, rule):
        if type(rule) == list:
            self.ruleStackruleStack = self.ruleStackruleStack + rule
        # TODO: Use error handling to check values passed, check if RuleFlood
        # TODO: Need to check for duplicate rules
        # TODO: Maybe use a set instead
        else:
            self.ruleStackruleStack.append(rule)
            
    def calc_delta(self, *, ruleType=RuleComplianceBase.MAX, tsOffset=0):
        # TODO: This is just a very simple implementation for Prado that just checks the defficit to max storage
        # TODO: This could be much more advanced to check diversions, controlled and uncontrolled releases
        # TODO: Also maxStor could be a guide curve rule or a scalar and have a check here which type
        if self.storMaxstorMax:
            storPrev = self.storstor[max(0,self.TT.step-1)].item()
            qInNet = self._get_qin_get_qin(tsOffset=tsOffset)
            lossEvap, storCur, qInNet = self._calc_evap_calc_evap(storPrev, qInNet)
            lossSeep, storCur, qInNet = self._calc_seepage_calc_seepage(storPrev, qInNet)
            qDiv, storCur, qInNet = self._calc_diversions_calc_diversions(storPrev, qInNet)
            return -max(0.,self.storMaxstorMax - storCur)*self.constantsconstants.af2cfs
        else: 
            return np.nan
 
    def _calc_storage(self, storPrev, qInNet, qOut):
        stor = storPrev - qOut*self.constantsconstants.cfs2af + qInNet*self.constantsconstants.cfs2af
        if stor <= 0.:
            qOut += stor*self.constantsconstants.af2cfs
            if qOut < 0.: qOut = 0.
            stor = 0.
        return stor, qOut
    
    def _calc_evap(self, stor, qIn):
        if self.evapObjevapObj:
            lossEvap = self.evapObjevapObj.calc_evap(stor+qIn*self.constantsconstants.cfs2af)*self.constantsconstants.af2cfs 
            stor, lossEvap = self._calc_storage_calc_storage(stor, qIn, lossEvap)
            qIn -= lossEvap
        else:
            lossEvap = 0.
        return lossEvap, stor, qIn
    
    def _calc_seepage(self, stor, qIn):
        if self.seepageObjseepageObj:
            lossSeep = self.seepageObjseepageObj.calc_seepage(stor+qIn*self.constantsconstants.cfs2af)
            stor, lossSeep = self._calc_storage_calc_storage(stor, qIn, lossSeep)
            qIn -= lossSeep
        else:
            lossSeep = 0.
        return lossSeep, stor, qIn
    
    def _calc_diversions(self, stor, qIn):
        qDiv = 0.
        if self.ruleDiversionruleDiversion:
            for curRule in self.ruleDiversionruleDiversion:
                curDiv = curRule.calc_release(rlsProposed=qIn+stor*self.constantsconstants.af2cfs,
                                              rlsPrev=curRule.release[max(0, self.TT.step-1)].item(),
                                              stor=stor, 
                                              qIn=qIn)
                storCur, curDiv, qIn = self._calc_storage_calc_storage(stor, qIn, curDiv)
                qIn -= curDiv
                qDiv += curDiv
        return qDiv, stor, qIn
 
    def calc_qout(self):
        qOutTotCur = 0.
        qInNet = qInInit = super()._get_qin()
        storPrev = self.storstor[max(0, self.TT.step-1)].item()
        storCur = storPrev
        # Get evaporation
        lossEvap, storCur, qInNet = self._calc_evap_calc_evap(storPrev, qInNet)
        if self.evapObjevapObj: self.lossEvaplossEvap[self.TT.step] = lossEvap
        # Get seepage
        lossSeep, storCur, qInNet = self._calc_seepage_calc_seepage(storPrev, qInNet)
        if self.seepageObjseepageObj: self.lossSeeplossSeep[self.TT.step] = lossSeep
        # Calc diversions
        qDiv, storCur, qInNet = self._calc_diversions_calc_diversions(storPrev, qInNet)
        if self.ruleDiversionruleDiversion: self.qDivqDiv[self.TT.step] = qDiv
        # Get uncontrolled release
        rlsUnCtrl = 0.
        # Get emergency release
        rlsEmgc = self.ruleEmgcruleEmgc.calc_release(
            0., stor=storPrev, rlsUnCtrl=rlsUnCtrl, qIn=qInNet) if self.ruleEmgcruleEmgc is not None else 0.
        qOutTotCur += rlsEmgc
        storCur, qOutTotCur = self._calc_storage_calc_storage(storPrev, qInNet, qOutTotCur)
        # TODO: If your emergency release exceeds your max controlled release then no rule stack
        # TODO: Should look into ruleStack as not a list but an ordered list and the rules unique
        # TODO: Python linked list seems like a good choice for rule stack
        qOutTotStart = qOutTotCur
        if self.outletCtrloutletCtrl:
            while True:
                rlsStack = self.calc_release_ctrlcalc_release_ctrl(self.rlsCtrlrlsCtrl[max(0, self.TT.step-1)],
                                                  rlsEmgc+rlsUnCtrl, storPrev, storCur, qInNet)
                storCur, qOutTotCur = self._calc_storage_calc_storage(storPrev, qInNet, qOutTotStart + rlsStack)
                rlsStack = qOutTotCur - qOutTotStart
                rlsUnCtrlPrev = rlsUnCtrl
                rlsUnCtrl = self.outletUnCtrloutletUnCtrl.calc_release(storCur) if self.outletUnCtrloutletUnCtrl else 0.
                qOutTotCur += rlsUnCtrl
                storCur, qOutTotCur = self._calc_storage_calc_storage(storPrev, qInNet, qOutTotCur)
                # TODO: This should also check if the ruleStack rules are type RLS_TOT if not then this check is not necessary
                if abs(rlsUnCtrl - rlsUnCtrlPrev) < 1e-10:
                    break
        else:
            rlsUnCtrl = self.outletUnCtrloutletUnCtrl.calc_release(storCur) if self.outletUnCtrloutletUnCtrl else 0.
            qOutTotCur += rlsUnCtrl
            storCur, qOutTotCur = self._calc_storage_calc_storage(storPrev, qInNet, qOutTotCur)
        if self.outletUnCtrloutletUnCtrl: self.rlsUnCtrlrlsUnCtrl[self.TT.step] = rlsUnCtrl
        # TODO: Check to see if a numpy reference is returned or a float primitive
        if self.outletCtrloutletCtrl: self.rlsCtrlrlsCtrl[self.TT.step] = min(rlsStack, self.outletCtrloutletCtrl.calc_release(storCur))
        self.qOutqOut[self.TT.step] = qOutTotCur
        self.continuitycontinuity[self.TT.step] = \
            (storCur-self.storstor[max(self.TT.step-1, 0)])*self.constantsconstants.af2cfs + \
                qOutTotCur - qInInit + lossEvap + lossSeep + qDiv
        self.storstor[self.TT.step] = max(0.,storCur)
        return qOutTotCur
    
    # TODO: This could be added to a rulestack object and generalized for diversions and release
    # TODO: RuleStack could subclass RuleBase for ulimate flexibility
    def calc_release_ctrl(self, rlsPrev, rlsUnCtrl, storPrev, storCur, qIn):
        rlsStack = 0.
        rlsMin = 0.
        rlsMax = np.inf
        # TODO: Should look into ruleStack as not a list but an ordered list and the rules unique
        # TODO: Python linked list seems like a good choice for rule stack
        # TODO: You could determine the rulease setting rule by determining which rule in the stack
        # set the release.
        if self.ruleStackruleStack:
            for curRule in self.ruleStackruleStack:
                rlsStack = curRule.calc_release(rlsStack,
                                                rlsPrev=rlsPrev,
                                                rlsUnCtrl=rlsUnCtrl,
                                                stor=storPrev, qIn=qIn)
                if hasattr(curRule, 'ruleType'):
                    if curRule.ruleType == RuleComplianceBase.MIN:
                        if curRule.qCompCur > rlsMin: rlsMin = curRule.qCompCur
                    elif curRule.ruleType == RuleComplianceBase.MAX:
                        if curRule.qCompCur < rlsMax: rlsMax = curRule.qCompCur
                storCur, rlsAdj = self._calc_storage_calc_storage(storPrev, qIn, rlsUnCtrl+rlsStack)
                rlsStack = rlsAdj - rlsUnCtrl
                if storCur == 0.: rlsMax = rlsStack
            # curRule.release[self.T.step] = rlsStack
            self.rlsMinrlsMin[self.TT.step] = rlsMin
            self.rlsMaxrlsMax[self.TT.step] = rlsMax
        return rlsStack
 
 
# NOTES:     
# Possible logic multiple spillways (not done)
     # def calc_release_unctrl(self, storInit):
     #     rlsUnCntrl = np.zeros(len(self.outletUnCtrl))
     #     rlsUnCntrlPrev = np.ones(len(self.outletUnCtrl))
     #     while all(rlsUnCntrl != rlsUnCntrlPrev):
     #         rlsUnCntrlPrev = rlsUnCntrl
     #         storCur = storInit
     #         for i, ctrlOutlet in self.outletUnCtrl:
     #             rlsUnCrl(i) = ctrlOutlet.calc_release
     #             storCur =