diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py new file mode 100644 index 0000000..ef58460 --- /dev/null +++ b/src/lisflood/hydrological_modules/mct.py @@ -0,0 +1,341 @@ +import numpy as np + + +def MCTRouting_single( + V00, q10, q01, q00, ql, q0mm, Cm0, Dm0, dt, xpix, s0, Balv, ANalv, Nalv +): + """ + This function implements Muskingum-Cunge-Todini routing method for a single pixel. + References: + Todini, E. (2007). A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach. Hydrol. Earth Syst. Sci. + (Chapter 5) + Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 + (Appendix B) + + :param V00: channel storage volume at t (beginning of computation step) + :param q10: O(t) - outflow (x+dx) at time t + :param q01: I(t+dt) - inflow (x) at time t+dt + :param q00: I(t) - inflow (x) at time t + :param ql: lateral flow over time dt [m3/s] + :param q0mm: I - average inflow during step dt + :param Cm0: Courant number at time t + :param Dm0: Reynolds number at time t + :param dt: time interval step + :param xpix: channel length + :param s0: channel slope + :param Balv: channel bankfull width + :param ANalv: angle of the riverbed side [rad] + :param Nalv: channel Manning roughness coefficient + :return: + q11: Outflow (x+dx) at O(t+dt) + V11: channel storage volume at t+dt + Cm1: Courant number at t+1 for state file + Dm1: Reynolds number at t+1 for state file + """ + + eps = 1e-06 + + # Calc O' first guess for the outflow at time t+dt + # O'(t+dt)=O(t)+(I(t+dt)-I(t)) + q11 = q10 + (q01 - q00) + + # check for negative and zero discharge values + # zero outflow is not allowed + if q11 < eps: # cmcheck <=0 + q11 = eps + + # calc reference discharge at time t + # qm0 = (I(t)+O(t))/2 + # qm0 = (q00 + q10) / 2. + + # Calc O(t+dt)=q11 at time t+dt using MCT equations + for i in range(2): # repeat 2 times for accuracy + + # reference I discharge at x=0 + qmx0 = (q00 + q01) / 2.0 + if qmx0 < eps: # cmcheck ==0 + qmx0 = eps + hmx0 = hoq(qmx0, s0, Balv, ANalv, Nalv) + + # reference O discharge at x=1 + qmx1 = (q10 + q11) / 2.0 + if qmx1 < eps: # cmcheck ==0 + qmx1 = eps + hmx1 = hoq(qmx1, s0, Balv, ANalv, Nalv) + + # Calc riverbed slope correction factor + cor = 1 - (1 / s0 * (hmx1 - hmx0) / xpix) + sfx = s0 * cor + if sfx < (0.8 * s0): + sfx = 0.8 * s0 # In case of instability raise from 0.5 to 0.8 + + # Calc reference discharge time t+dt + # Q(t+dt)=(I(t+dt)+O'(t+dt))/2 + qm1 = (q01 + q11) / 2.0 + # cm + if qm1 < eps: # cmcheck ==0 + qm1 = eps + # cm + hm1 = hoq(qm1, s0, Balv, ANalv, Nalv) + dummy, Ax1, Bx1, Px1, ck1 = qoh(hm1, s0, Balv, ANalv, Nalv) + if ck1 <= eps: + ck1 = eps + + # Calc correcting factor Beta at time t+dt + Beta1 = ck1 / (qm1 / Ax1) + # calc corrected cell Reynolds number at time t+dt + Dm1 = qm1 / (sfx * ck1 * Bx1 * xpix) / Beta1 + # corrected Courant number at time t+dt + Cm1 = ck1 * dt / xpix / Beta1 + + # Calc MCT parameters + den = 1 + Cm1 + Dm1 + c1 = (-1 + Cm1 + Dm1) / den + c2 = (1 + Cm0 - Dm0) / den * (Cm1 / Cm0) + c3 = (1 - Cm0 + Dm0) / den * (Cm1 / Cm0) + c4 = (2 * Cm1) / den + + # cmcheck + # Calc outflow q11 at time t+1 + # Mass balance equation without lateral flow + # q11 = c1 * q01 + c2 * q00 + c3 * q10 + # Mass balance equation that takes into consideration the lateral flow + q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql + + if q11 < eps: # cmcheck <=0 + q11 = eps + + #### end of for loop + + # # cmcheck + # calc_t = xpix / ck1 + # if calc_t < dt: + # print('xpix/ck1 < dt') + + # k1 = dt / Cm1 + # x1 = (1. - Dm1) / 2. + + # Calc the corrected mass-conservative expression for the reach segment storage at time t+dt + # The lateral inflow ql is only explicitly accounted for in the mass balance equation, while it is not in the equation expressing + # the storage as a weighted average of inflow and outflow.The rationale of this approach lies in the fact that the outflow + # of the reach implicitly takes the effect of the lateral inflow into account. + if q11 == 0: + V11 = V00 + (q00 + q01 - q10 - q11) * dt / 2 + else: + V11 = (1 - Dm1) * dt / (2 * Cm1) * q01 + (1 + Dm1) * dt / (2 * Cm1) * q11 + # V11 = k1 * (x1 * q01 + (1. - x1) * q11) # MUST be the same as above! + + ### calc integration on the control volume (pixel) + # calc average discharge outflow q1m for MCT channels during routing sub step dt + # Calculate average outflow using water balance for MCT channel grid cell over sub-routing step + q1mm = q0mm + ql + (V00 - V11) / dt + + # cmcheck + # q1m cannot be smaller than eps or it will cause instability + if q1mm < eps: + q1mm = eps + V11 = V00 + (q0mm + ql - q1mm) * dt + + # q11 Outflow at O(t+dt) + # q1m average outflow in time dt + # V11 water volume at t+dt + # Cm1 Courant number at t+dt for state files + # Dm1 Reynolds numbers at t+dt for state files + return q11, V11, q1mm, Cm1, Dm1 + + +def hoq(q, s0, Balv, ANalv, Nalv): + """Water depth h from discharge q. + Given a generic cross-section (rectangular, triangular or trapezoidal) and a steady-state discharge q=Q*, it computes + water depth (y), wet contour (Bx), wet area (Ax) and wave celerity (cel) using Newton-Raphson method. + Reference: + Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. + Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 + + :param q: steady-state river discharge [m3/s] + :param s0: river bed slope (tan B) + :param Balv : width of the riverbed [m] + :param ChanSdXdY : slope dx/dy of riverbed side + :param ANalv : angle of the riverbed side [rad] + :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] + + :return: + y: water depth referred to the bottom of the riverbed [m] + """ + + alpha = 5.0 / 3.0 # exponent (5/3) + eps = 1.0e-06 + max_tries = 1000 + + rs0 = np.sqrt(s0) + usalpha = 1.0 / alpha + + # cotangent(angle of the riverbed side - dXdY) + if ANalv < np.pi / 2: + # triangular or trapezoid cross-section + c = cotan(ANalv) + else: + # rectangular corss-section + c = 0.0 + + # sin(angle of the riverbed side - dXdY) + if ANalv < np.pi / 2: + # triangular or trapezoid cross-section + s = np.sin(ANalv) + else: + # rectangular cross-section + s = 1.0 + + # water depth first approximation y0 based on steady state q + if Balv == 0: + # triangular cross-section + y = (Nalv * q / rs0) ** (3.0 / 8.0) * (2 / s) ** 0.25 / c ** (5.0 / 8.0) + else: + # rectangular cross-section and first approx for trapezoidal cross-section + y = (Nalv * q / (rs0 * Balv)) ** usalpha + + if (Balv != 0) and (ANalv < np.pi / 2): + # trapezoid cross-section + y = (Nalv * q / rs0) ** usalpha * (Balv + 2.0 * y / s) ** 0.4 / (Balv + c * y) + + for tries in range(1, max_tries): + # calc Q(y) for the different tries of y + q0, Ax, Bx, Px, cel = qoh(y, s0, Balv, ANalv, Nalv) + # Ax: wet area[m2] + # Bx: cross-section width at water surface[m] + # Px: cross-section wet contour [m] + # cel: wave celerity[m/s] + + # this is the function we want to find the 0 for f(y)=Q(y)-Q* + fy = q0 - q + # calc first derivative of f(y) f'(y)=Bx(y)*cel(y) + dfy = Bx * cel + # calc update for water depth y + dy = fy / dfy + # update yt+1=yt-f'(yt)/f(yt) + y = y - dy + # stop loop if correction becomes too small + if np.abs(dy) < eps: + break + + return y + + +def qoh(y, s0, Balv, ANalv, Nalv): + """Discharge q from water depth h. + Given a generic river cross-section (rectangular, triangular and trapezoidal) + and a water depth (y [m]) referred to the bottom of the riverbed, it uses Manning’s formula to calculate: + q: steady-state discharge river discharge [m3/s] + a: wet area [m2] + b: cross-section width at water surface [m] + p: cross-section wet contour [m] + cel: wave celerity [m/s] + Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 + + :param y: river water depth [m] + :param s0: river bed slope (tan B) + :param Balv : width of the riverbed [m] + :param ANalv : angle of the riverbed side [rad] + :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] + + :return: + q,a,b,p,cel + """ + + alpha = 5.0 / 3.0 # exponent (5/3) + + rs0 = np.sqrt(s0) + alpham = alpha - 1.0 + + # np.where(ANalv < np.pi/2, triang. or trapeiz., rectangular) + # cotangent(angle of the riverbed side - dXdY) + c = np.where( + ANalv < np.pi / 2, + # triangular or trapezoid cross-section + cotan(ANalv), + # rectangular cross-section + 0.0, + ) + # sin(angle of the riverbed side - dXdY) + s = np.where( + ANalv < np.pi / 2, + # triangular or trapezoid cross-section + np.sin(ANalv), + # rectangular corss-section + 1.0, + ) + + a = (Balv + y * c) * y # wet area [m2] + b = Balv + 2.0 * y * c # cross-section width at water surface [m] + p = Balv + 2.0 * y / s # cross-section wet contour [m] + q = rs0 / Nalv * a**alpha / p**alpham # steady-state discharge [m3/s] + cel = (q / 3.0) * (5.0 / a - 4.0 / (p * b * s)) # wave celerity [m/s] + + return q, a, b, p, cel + + +def hoV(V, xpix, Balv, ANalv): + """Water depth h from volume V. + Given a generic river cross-section (rectangular, triangular and trapezoidal) and a river channel volume V, + it calculates the water depth referred to the bottom of the riverbed [m] (y). + Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 + + :param V : volume of water in channel riverbed + :param xpix : dimension along the flow direction [m] + :param Balv : width of the riverbed [m] + :param ANalv : angle of the riverbed side [rad] + + :return: + y : channel water depth [m] + """ + eps = 1e-6 + + c = np.where( + ANalv < np.pi / 2, # angle of the riverbed side dXdY [rad] + cotan(ANalv), # triangular or trapezoidal cross-section + 0.0, + ) # rectangular cross-section + + a = V / xpix # wet area [m2] + + # np.where(c < 1.d-6, rectangular, triangular or trapezoidal) + y = np.where( + np.abs(c) < eps, + a / Balv, # rectangular cross-section + (-Balv + np.sqrt(Balv**2 + 4 * a * c)) / (2 * c), + ) # triangular or trapezoidal cross-section + + return y + + +def qoV(V, xpix, s0, Balv, ANalv, Nalv): + """Discharge q from river channel volume V. + Given a generic river cross-section (rectangular, triangular and trapezoidal) + and a water volume (V [m3]), it uses Manning’s formula to calculate the corresponding discharge (q [m3/s]). + Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 + + :param V : volume of water in channel riverbed + :param xpix : dimension along the flow direction [m] + :param s0: river bed slope (tan B) + :param Balv : width of the riverbed [m] + :param ANalv : angle of the riverbed side [rad] + :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] + + :return: + y : channel water depth [m] + """ + y = hoV(V, xpix, Balv, ANalv) + q, a, b, p, cel = qoh(y, s0, Balv, ANalv, Nalv) + return q + + +def cotan(x): + """There is no cotangent function in numpy""" + return np.cos(x) / np.sin(x) + + +def rad_from_dxdy(dxdy): + """Calculate radians""" + rad = np.arctan(1 / dxdy) + # angle = np.rad2deg(rad) + return rad diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index c3c4b5c..50bbb62 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -27,6 +27,7 @@ from .inflow import inflow from .transmission import transmission from .kinematic_wave_parallel import kinematicWave, kwpt +from .mct import MCTRouting_single from ..global_modules.settings import LisSettings, MaskInfo from ..global_modules.add1 import loadmap, loadmap_base, compressArray, decompress @@ -376,7 +377,7 @@ def initial(self): # cmcheck # initialising the variable # should initialisation be moved somewhere else? - self.var.ChanQAvgDt = maskinfo.in_zero() + self.var.ChanQAvgDt = maskinfo.in_zero() self.var.ChanQKinAvgDt = maskinfo.in_zero() # cmcheck @@ -404,16 +405,13 @@ def initial(self): # cmcheck - non so se sostituita da self.var.sumInWB self.var.sumInWB = maskinfo.in_zero() - def initialSecond(self): """ initial part of the second line channel routing module for split routing """ settings = LisSettings.instance() option = settings.options - binding = settings.binding flags = settings.flags - self.var.ChannelAlpha2 = None # Default value, if split-routing is not active and water is routed only in the main channel @@ -523,7 +521,6 @@ def initialMCT(self): """ settings = LisSettings.instance() option = settings.options - binding = settings.binding if not (option['SplitRouting']): self.var.ChannelAlpha2 = None @@ -548,7 +545,6 @@ def initialMCT(self): # Identify channel pixels where Kinematic wave is used instead of MCT self.var.LddMCT = lddmask(self.var.LddChan, self.var.IsChannelMCTPcr) #pcr - self.var.LddMCTNp = compressArray(self.var.LddMCT) #np # Ldd for MCT routing self.var.LddKinematic = lddmask(self.var.LddChan, self.var.IsChannelKinematicPcr) #pcr @@ -594,7 +590,7 @@ def initialMCT(self): # ************************************************************ # ***** INITIALISE MUSKINGUM-CUNGE-TODINI WAVE ROUTER ******** # ************************************************************ - self.mct_river_router = mctWave(self.get_mct_pix(compressArray(self.var.LddMCT)), self.var.mctmask) + self.mct_river_router = mctWave(self.compress_mct(compressArray(self.var.LddMCT)), self.var.mctmask) # -------------------------------------------------------------------------- # -------------------------------------------------------------------------- @@ -604,7 +600,6 @@ def dynamic(self, NoRoutingExecuted): """ settings = LisSettings.instance() option = settings.options - binding = settings.binding if not(option['InitLisflood']): # only with no InitLisflood self.lakes_module.dynamic_inloop(NoRoutingExecuted) @@ -675,26 +670,20 @@ def dynamic(self, NoRoutingExecuted): if option['wateruse']: self.var.AddedTRUN -= np.take(np.bincount(self.var.Catchments, weights=self.var.WUseAddM3Dt.copy()),self.var.Catchments) - - SideflowChan = np.where(self.var.IsChannelKinematic, SideflowChanM3 * self.var.InvChanLength * self.var.InvDtRouting,0) # Sideflow contribution to kinematic and split routing grid cells expressed in [cu m /s / m channel length] - - #cmcheck - # non c'e' bisogno di questo if? - if option['MCTRouting']: - SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting,0) #Ql - # Sideflow contribution to MCT grid cells expressed in [m3/s] - else: - SideflowChanMCT = 0 + SideflowChan = np.where(self.var.IsChannelKinematic, SideflowChanM3 * self.var.InvChanLength * self.var.InvDtRouting,0) # ************************************************************ # ***** RIVER ROUTING **************** # ************************************************************ - if option['InitLisflood']: self.var.IsChannelKinematic = self.var.IsChannel.copy() + if option['InitLisflood']: + self.var.IsChannelKinematic = self.var.IsChannel.copy() # Use kinematic routing in all grid cells - # KINEMATIC ROUTING ONLY - InitLisflood - if option['InitLisflood'] or (not(option['SplitRouting']) and (not(option['MCTRouting']))): + # only run kinematic for InitLisflood option + + # KINEMATIC ROUTING ONLY + if option['InitLisflood'] or not (option['SplitRouting']): # Kinematic routing uses the same variables for inputs and outputs. Values are updated during computation # from upstream to downstream. @@ -702,25 +691,17 @@ def dynamic(self, NoRoutingExecuted): # Store results to general variables for channels routing - self.var.ChanQ = self.var.ChanQKin.copy() + ChanQ = self.var.ChanQKin.copy() # Outflow (x+dx) Q at the end of computation step t+dt for full section (instant) - self.var.ChanM3 = self.var.ChanM3Kin.copy() + ChanM3 = self.var.ChanM3Kin.copy() # Channel storage V at the end of computation step t+dt for full section (instant) - self.var.ChanQAvgDt = self.var.ChanQKinAvgDt.copy() + ChanQAvgDt = self.var.ChanQKinAvgDt.copy() # Average Outflow (x+dx) Q during computation step dt for full section (instant) - # self.var.sumDisDay += self.var.ChanQ - # sum of river outflow on model sub-steps (instantaneous discharge) - - self.var.sumDisDay += self.var.ChanQAvgDt - # sum of river outflow on model sub-steps (average discharge) - # This is used to calculate the average outflow discharge on the model time step - - - # SPLIT ROUTING ONLY - no InitLisflood - if not option['InitLisflood'] and option['SplitRouting'] and not(option['MCTRouting']): + # SPLIT ROUTING ONLY + else: # Split routing uses the same variables for inputs and outputs. Values are updated during computation # from upstream to downstream. @@ -729,136 +710,84 @@ def dynamic(self, NoRoutingExecuted): # Combine the two lines of routing together # Store results to general variables for channels routing - self.var.ChanQ = np.maximum(self.var.ChanQKin + self.var.Chan2QKin - self.var.QLimit, 0) + ChanQ = np.maximum(self.var.ChanQKin + self.var.Chan2QKin - self.var.QLimit, 0) # (real) total outflow (at x+dx) at time t+dt end of step for the full cross-section (instant) # Main channel routing and above bankfull routing from second line of routing - self.var.ChanM3 = self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start + ChanM3 = self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start # Total channel storage [m3] = Volume in main channel (ChanM3Kin) + volume above bankfull in second line (Chan2M3Kin - Chan2M3Start) # Total channel storage V at the end of computation step t+dt for full section (instant) # cmcheck - self.var.ChanQAvgDt = np.maximum(self.var.ChanQKinAvgDt + self.var.Chan2QKinAvgDt - self.var.QLimit, 0) + ChanQAvgDt = np.maximum(self.var.ChanQKinAvgDt + self.var.Chan2QKinAvgDt - self.var.QLimit, 0) # (real) total outflow (at x+dx) at time t+dt end of step for the full cross-section (instant) # Main channel routing and above bankfull routing from second line of routing - # cmcheck - # self.var.sumDisDay += self.var.ChanQ - # sum of total river outflow on model sub-step (instantaneous discharge) - - self.var.sumDisDay += self.var.ChanQAvgDt - # sum of total river outflow on model sub-step (instantaneous discharge) - # This is used to calculate the average outflow discharge on the model time step - - # KINEMATIC/SPLIT ROUTING AND MUSKINGUM-CUNGE-TODINI - no InitLisflood - if not option['InitLisflood'] and option['MCTRouting']: + if option['MCTRouting']: # To parallelise the routing, pixels are grouped by orders. Pixels in the same order are independent and # routing can be solved in parallel. # The same order can have both Kin and MCT pixels. # First, Kinematic routing is solved on all pixels (including MCT pixels) to generate the inputs to the # MCT grid cells downstream (MCT needs input from upstream kinematic pixels). # (this needs to be changed in future versions) - - if not (option['SplitRouting']): - # Kinematic routing (+ MCT) - # This is calculated for every grid cell, including MCT grid cells - - ### calling kinematic routing - self.KinRouting(SideflowChan) - - # saving outputs to be used as inputs for MCT routing - # Do not use self.var.ChanQ and self.var.ChanM3 because they store the output of MCT and it must - # not be over-written by the output of the kinematic routing. - - ChanQKinOutEnd = self.var.ChanQKin.copy() - # Outflow (x+dx) Q at the end of computation step t+dt for full section (instant) - - ChanM3KinEnd = self.var.ChanM3Kin.copy() - # Channel storage V at the end of computation step t+dt for full section (instant) - - ChanQKinOutAvgDtEnd = self.var.ChanQKinAvgDt.copy() - # Average Outflow (x+dx) Q during computation step dt for full section (instant) - - - if (option['SplitRouting']): - # Split routing (+ MCT) - # This is calculated for every grid cell, including MCT grid cells - - ### calling split routing - self.SplitRouting(SideflowChan) - - # Combine the two lines of routing together - # Saving outputs to be used as inputs for MCT routing - # Do not use self.var.ChanQ and self.var.ChanM3 because they store the output of MCT and it must - # not be over-written by the output of the split routing. - - ChanQKinOutEnd = np.maximum(self.var.ChanQKin + self.var.Chan2QKin - self.var.QLimit, 0) - # (real) Total outflow (at x + dx) at time t + dt (instant) - # Main channel routing and above bankfull routing from second line of routing - - ChanM3KinEnd = self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start - # (real) Total channel storage (at x + dx) at time t + dt (instant) - # Total channel storage [m3] = Volume in main channel (ChanM3Kin) + volume above bankfull in second - # line (Chan2M3Kin - Chan2M3Start) - - # cmcheck - ChanQKinOutAvgDtEnd = np.maximum(self.var.ChanQKinAvgDt + self.var.Chan2QKinAvgDt - self.var.QLimit,0) - # (real) total outflow (at x+dx) at time t+dt end of step for the full cross-section (instant) - # Main channel routing and above bankfull routing from second line of routing - + + # Sideflow contribution to MCT grid cells expressed in [m3/s] + SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql # MCT routing # This is calculated for MCT grid cell only but takes the output of kinematic or split routing - ChanQMCTOutStart = self.var.ChanQ.copy() + # ChanQMCTOutStart = self.var.ChanQ.copy() # Outflow (x+dx) at time t (instant) q10 # cmcheck - debug # ChanQKinOutEnd[ChanQKinOutEnd != 0] = 0 # Set inflow from kinematic pixels to 1 - ChanM3MCTStart = self.var.ChanM3.copy() + # ChanM3MCTStart = self.var.ChanM3.copy() # Channel storage at time t V00 - ChanQMCTInStart = self.var.PrevQMCTin.copy() + # ChanQMCTInStart = self.var.PrevQMCTin.copy() # Inflow (x) at time t instant (instant) q00 # This is coming from upstream pixels #### calling MCT routing # Using ChanQKinOutEnd and ChanQKinOutAvgDtEnd as inputs to receive inflow from upstream kinematic/split # routing pixels that are contributing to MCT pixels - ChanQMCTOutAvgDt,ChanQMCTOutEnd,ChanM3MCTEnd,Cmend, Dmend = self.MCTRoutingLoop(ChanQMCTOutStart,ChanQMCTInStart,ChanQKinOutEnd,ChanQKinOutAvgDtEnd,SideflowChanMCT,ChanM3MCTStart) - # Average outflow (x+dx) during time dt step (average) q1m - # Outflow (x+dx) at time t+dt end of calculation step (instant) q11 - # Channel storage at time t+dt end of calculation step (instant) V11 - - # Storing MCT Courant and Reynolds numbers for state files - self.var.PrevCm0 = Cmend - self.var.PrevDm0 = Dmend - - # combine results from Kinematic and MCT pixels - # Store results to general variables for channels routing - - self.var.ChanQ = np.where(self.var.IsChannelKinematic, ChanQKinOutEnd, ChanQMCTOutEnd) - # Outflow (x+dx) Q at the end of computation step t+dt (instant) - self.var.ChanM3 = np.where(self.var.IsChannelKinematic, ChanM3KinEnd, ChanM3MCTEnd) - # Channel storage V at the end of computation step t+dt (instant) - self.var.ChanQAvgDt = np.where(self.var.IsChannelKinematic, ChanQKinOutAvgDtEnd, ChanQMCTOutAvgDt) - # Average channel outflow (x+dx) QAvg during routing step (average) - - # update input (x) Q at t for next step (instant) - ChanQMCTStartPcr = decompress(self.var.ChanQ) # pcr - self.var.PrevQMCTin = compressArray(upstream(self.var.LddChan, ChanQMCTStartPcr)) - # using LddChan here because we need to include the input from upstream kinematic pixels - - self.var.sumDisDay += self.var.ChanQAvgDt - # sum of average outflow Q on model sub-steps to give accumulated total outflow volume - - # This should NOT be used - # self.var.sumDisDay += self.var.ChanQ - # sum of total river outflow on model sub-step + + # Grab current state + ChanQ_0 = self.var.ChanQ.copy() + ChanQAvgDt_0 = self.var.ChanQAvgDt.copy() + ChanM3_0 = self.var.ChanM3.copy() + + # Grab next state, but incomplete (only Kinematic/Split points) + ChanQ_1 = ChanQ + ChanQAvgDt_1 = ChanQAvgDt + + # put results of Kinematic/Split routing into MCT points + self.var.ChanQ = ChanQ + self.var.ChanM3 = ChanM3 + self.var.ChanQAvgDt = ChanQAvgDt + + # Update current state at MCT points + self.MCTRoutingLoop( + ChanQ_0, + ChanQAvgDt_0, + ChanM3_0, + ChanQ_1, + ChanQAvgDt_1, + SideflowChanMCT + ) + else: + # put results of Kinematic/Split routing into MCT points + self.var.ChanQ = ChanQ + self.var.ChanM3 = ChanM3 + self.var.ChanQAvgDt = ChanQAvgDt + + # sum of average outflow Q on model sub-steps to give accumulated total outflow volume + self.var.sumDisDay += self.var.ChanQAvgDt ### end of river routing calculation # ---- Uncomment lines 603-635 in order to compute the mass balance error within the routing module for the options (i) initial run or (ii) split routing off ---- @@ -1129,18 +1058,27 @@ def SplitRouting(self, SideflowChan): # (virtual) total outflow (at x + dx) at time t + dt (instant) for second line of routing (using increased Manninig coeff) # Correct negative discharge at the end of computation step in second line - FldpQKin = self.var.Chan2QKin - self.var.QLimit + # FldpQKin = self.var.Chan2QKin - self.var.QLimit # Outflow at t+dt from above bankfull only - return + def upstream(self, var): + var_pcr = decompress(var) + var_ups = compressArray(upstream(self.var.LddChan, var_pcr)) + return var_ups - def MCTRoutingLoop(self,ChanQMCTOutStart,ChanQMCTInStart,ChanQKinOut,ChanQKinOutAvgDt,SideflowChanMCT,ChanM3Start): + def MCTRoutingLoop( + self, + ChanQ_0, + ChanQAvgDt_0, + ChanM3_0, + ChanQ_1, + SideflowChanMCT + ): """This function implements Muskingum-Cunge-Todini routing method MCT routing is calculated on MCT pixels only but gets inflow from both KIN and MCT upstream pixels. This is the reason we need both ChanQKinOutEnd and ChanQKinOutAvgDt as inputs. - Function get_mct_pix is used to compress arrays with all river channel pixels to arrays containing MCT pixels only. - Function put_mct_pix is used to explode arrays containing MCT pixels only back to arrays containing all rivers pixels. + Function compress_mct is used to compress arrays with all river channel pixels to arrays containing MCT pixels only. References: Todini, E. (2007). A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach. Hydrol. Earth Syst. Sci. (Chapter 5) @@ -1161,83 +1099,13 @@ def MCTRoutingLoop(self,ChanQMCTOutStart,ChanQMCTInStart,ChanQKinOut,ChanQKinOut Dmout: Reynolds number at time t+dt """ - # channel geometry - MCT pixels only - xpix = self.get_mct_pix(self.var.ChanLength) # dimension along the flow direction [m] - s0 = self.get_mct_pix(self.var.ChanGrad) # river bed slope (tan B) - Balv = self.get_mct_pix(self.var.ChanBottomWidth) # width of the riverbed [m] - ChanSdXdY = self.get_mct_pix(self.var.ChanSdXdY) # slope dx/dy of riverbed side - # Nalv = self.get_mct_pix(self.var.ChanMan) # channel mannings coefficient n for the riverbed [s/m1/3] - Nalv = self.get_mct_pix(self.var.ChanManMCT) # channel mannings coefficient n for the riverbed [s/m1/3] - ANalv = np.arctan(1/ChanSdXdY) # angle of the riverbed side [rad] - - dt = self.var.DtSecChannel # computation time step for channel [s] - - # MCT Courant and Reynolds numbers from previous step (MCT pixels) - Cm0 = self.get_mct_pix(self.var.PrevCm0) - Dm0 = self.get_mct_pix(self.var.PrevDm0) - - # instant discharge at channel input (I x=0) and channel output (O x=1) - # at the end of previous calculation step I(t) and O(t) and - # at the end of current calculation step I(t+1) and O(t+1) - - # Inflow at time t - # I(t) - # calc contribution from upstream pixels at time t (instant) (dim=all pixels) - q00 = self.get_mct_pix(ChanQMCTInStart) - - V00 = self.get_mct_pix(ChanM3Start) - # Water volume in the channel at beginning of computation step (t) (and end of previous time step) - - # Outflow at time t - # O(t) - # dim=mct pixels - q10 = self.get_mct_pix(ChanQMCTOutStart) - - # calc contribution from upstream pixels at time t+1 (dim=all pixels because we need to include both MCT and KIN pixels at the same time) - # First order (aka order 0) of mct pixels have inflow and average inflow from kin pixels only - ChanQMCTPcr=decompress(ChanQKinOut) #pcr - ChanQMCTUp1=compressArray(upstream(self.var.LddChan,ChanQMCTPcr)) - # Inflow at time t+1 - # I(t+dt) - # dim=mct pixels - q01 = self.get_mct_pix(ChanQMCTUp1) - - # calc average discharge from upstream pixels during step dt - # First order (aka order 0) of mct pixels have inflow and average inflow from kin pixels only - # I have already calculated outflow for all kin pixels for this step - ChanQKinOutAvgDtPcr=decompress(ChanQKinOutAvgDt) #pcr - ChanQKinOutAvgDtUp1=compressArray(upstream(self.var.LddChan,ChanQKinOutAvgDtPcr)) - # Average Inflow during time dt - # I average (dt) - # dim=mct pixels - q0m = self.get_mct_pix(ChanQKinOutAvgDtUp1) - - # Outflow (x+1) at time t+1 - # O(t+dt) - # dim=mct pixels - # set to zero at beginning of computation - q11 = np.zeros_like(q01) - q1m = np.zeros_like(q01) - V11 = np.zeros_like(q01) - - # Lateral flow Ql (average) during interval dt [m3/s] - # Ql(t) - # calc contribution from sideflow - ql = self.get_mct_pix(SideflowChanMCT) - - - ### start pixels loop on all mct pixels ### - # Pixels in the same order are independent and can be routed in parallel. - # Orders must be processed in series, starting from order 0. - # Outflow from MCT pixels can only go to a MCT pixel - # First order (aka order 0) of mct pixels have inflow and average inflow from kin pixels only - ChanQOut = ChanQKinOut.copy() - ChanQOutAvg = ChanQKinOutAvgDt.copy() - # Initialise outflow (x+dx) at t+dt using outflow from kinematic pixels - # Then update outflow from mct pixels as different orders are calculated - # Same for the average outflow + dt = self.var.DtSecChannel # computation time step for channel [s] num_orders = self.mct_river_router.order_start_stop.shape[0] + + index_kin = range(len(self.var.ChanQ)) + mapping_mct = self.compress_mct(index_kin) + # loop on orders for order in range(num_orders): first = self.mct_river_router.order_start_stop[order, 0] @@ -1245,400 +1113,47 @@ def MCTRoutingLoop(self,ChanQMCTOutStart,ChanQMCTInStart,ChanQKinOut,ChanQKinOut # loop on pixels in the order for index in range(first, last): # get pixel ID - idpix = self.mct_river_router.pixels_ordered[index] - - # cmcheck - eps = 1e-06 - # check upstream inflow - # q01 cannot be zero - # This is a problem for head pixels - # if q01[idpix] == 0: - # q01[idpix] = eps - - ### calling MCT function for single cell idpx - q11[idpix], V11[idpix],q1m[idpix], Cm0[idpix], Dm0[idpix] = self.MCTRouting_single(V00[idpix],q10[idpix], q01[idpix], q00[idpix], ql[idpix],q0m[idpix], Cm0[idpix], Dm0[idpix], - dt, xpix[idpix], s0[idpix], Balv[idpix], ANalv[idpix], Nalv[idpix]) - # # debug - cmcheck - tanto entra tanto esce - # q11[idpix] = q01[idpix] + ql[idpix] - # # Set outflow to be the same as inflow in MCT grid cells ('tanto entra tanto esce') - # V11[idpix] = V00[idpix] - # # Set end volume to be the same as initial volume in MCT grid cells ('tanto entra tanto esce') - - # Update contribution from upstream pixels at time t+1 (dim=all pixels) using the newly calculated q11 - # I want to update q01 (inflow at t+1) for cells downstream of idpix using the newly calculated q11 - # It can be a mix of kinematic and mct pixels. I want to update the contribution from mct pixels - Q11 = self.put_mct_pix(q11) - # explode mct pixels to all catchment pixels - - # Update average contribution from upstream pixels at time t+1 (dim=all pixels) using the newly calculated q1m - # I want to update q0m (avergae inflow at t+1) for cells downstream of this order using the newly calculated q1m, because - # downstream pixels receive outflow from this order and from mct pixels in this order. - # It can be a mix of kinematic and mct pixels. I want to update the contribution from mct pixels - Q1m = self.put_mct_pix(q1m) - # explode mct pixels to all catchment pixels - - # combine results in pixels from this order (mct pixels) with results in pixels from upstream orders (kin and mct) - # pixels that are not MCT get value zero in Q11 from put_mct_pix (see above) - ChanQOut = np.where(Q11 == 0, ChanQOut, Q11) - ChanQOutAvg = np.where(Q1m == 0, ChanQOutAvg, Q1m) - - # for each pixel in the catchment, calc contribution from upstream pixels - # outflow (x+dx) at t+dt for both kin and mct upstream pixels is stored in ChanQOut - QupPcr=decompress(ChanQOut) #pcr - Qup01=compressArray(upstream(self.var.LddChan,QupPcr)) - - # for each pixel in the catchment, calc average contribution from upstream pixels - # average outflow (x+dx) at t+dt for both kin and mct upstream pixels is stored in ChanQOutAvg - QupAvgPcr=decompress(ChanQOutAvg) #pcr - Qup0m=compressArray(upstream(self.var.LddChan,QupAvgPcr)) - - # slice the MCT pixels - # Inflow at time t+1 - # I(t+dt) - # dim=mct pixels - q01 = self.get_mct_pix(Qup01) - q0m = self.get_mct_pix(Qup0m) - - # repeat for next order of pixels - - ### end order loop - all pixels have been calculated ### - - # explode arrays with MCT results on all catchment pixels - ChanQMCTOut = self.put_mct_pix(q11) - Cmout = self.put_mct_pix(Cm0) - Dmout = self.put_mct_pix(Dm0) - ChanM3MCTOut = self.put_mct_pix(V11) - ChanQMCTOutAvg = self.put_mct_pix(q1m) - - return ChanQMCTOutAvg, ChanQMCTOut, ChanM3MCTOut, Cmout, Dmout - - - def MCTRouting_single(self, V00, q10, q01, q00, ql,q0mm, Cm0, Dm0, dt, xpix, s0, Balv, ANalv, Nalv): - """ - This function implements Muskingum-Cunge-Todini routing method for a single pixel. - References: - Todini, E. (2007). A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach. Hydrol. Earth Syst. Sci. - (Chapter 5) - Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 - (Appendix B) - - :param V00: channel storage volume at t (beginning of computation step) - :param q10: O(t) - outflow (x+dx) at time t - :param q01: I(t+dt) - inflow (x) at time t+dt - :param q00: I(t) - inflow (x) at time t - :param ql: lateral flow over time dt [m3/s] - :param q0mm: I - average inflow during step dt - :param Cm0: Courant number at time t - :param Dm0: Reynolds number at time t - :param dt: time interval step - :param xpix: channel length - :param s0: channel slope - :param Balv: channel bankfull width - :param ANalv: angle of the riverbed side [rad] - :param Nalv: channel Manning roughness coefficient - :return: - q11: Outflow (x+dx) at O(t+dt) - V11: channel storage volume at t+dt - Cm1: Courant number at t+1 for state file - Dm1: Reynolds number at t+1 for state file - """ - - eps = 1e-06 - - # Calc O' first guess for the outflow at time t+dt - # O'(t+dt)=O(t)+(I(t+dt)-I(t)) - q11 = q10 + (q01 - q00) - - # check for negative and zero discharge values - # zero outflow is not allowed - if q11 < eps: # cmcheck <=0 - q11 = eps - - # calc reference discharge at time t - # qm0 = (I(t)+O(t))/2 - # qm0 = (q00 + q10) / 2. - - # Calc O(t+dt)=q11 at time t+dt using MCT equations - for i in range(2): # repeat 2 times for accuracy - - # reference I discharge at x=0 - qmx0 = (q00 + q01) / 2. - if qmx0 < eps: # cmcheck ==0 - qmx0 = eps - hmx0 = self.hoq(qmx0, s0, Balv, ANalv, Nalv) - - # reference O discharge at x=1 - qmx1 = (q10 + q11) / 2. - if qmx1 < eps: # cmcheck ==0 - qmx1 = eps - hmx1 = self.hoq(qmx1, s0,Balv,ANalv,Nalv) - - # Calc riverbed slope correction factor - cor = 1 - (1 / s0 * (hmx1 - hmx0) / xpix) - sfx = s0 * cor - if sfx < (0.8 * s0): - sfx = 0.8 * s0 # In case of instability raise from 0.5 to 0.8 - - # Calc reference discharge time t+dt - # Q(t+dt)=(I(t+dt)+O'(t+dt))/2 - qm1 = (q01 + q11) / 2. - #cm - if qm1 < eps: # cmcheck ==0 - qm1 = eps - #cm - hm1 = self.hoq(qm1,s0,Balv,ANalv,Nalv) - dummy, Ax1,Bx1,Px1,ck1 = self.qoh(hm1,s0,Balv,ANalv,Nalv) - if (ck1 <= eps): - ck1 = eps - - # Calc correcting factor Beta at time t+dt - Beta1 = ck1 / (qm1 / Ax1) - # calc corrected cell Reynolds number at time t+dt - Dm1 = qm1 / (sfx * ck1 * Bx1 * xpix) / Beta1 - # corrected Courant number at time t+dt - Cm1 = ck1 * dt / xpix / Beta1 - - # Calc MCT parameters - den = 1 + Cm1 + Dm1 - c1 = (-1 + Cm1 + Dm1) / den - c2 = (1 + Cm0 - Dm0) / den * (Cm1 / Cm0) - c3 = (1 - Cm0 + Dm0) / den * (Cm1 / Cm0) - c4 = (2 * Cm1) / den - - # cmcheck - # Calc outflow q11 at time t+1 - # Mass balance equation without lateral flow - # q11 = c1 * q01 + c2 * q00 + c3 * q10 - # Mass balance equation that takes into consideration the lateral flow - q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql - - if q11 < eps: # cmcheck <=0 - q11 = eps - - #### end of for loop - - # # cmcheck - # calc_t = xpix / ck1 - # if calc_t < dt: - # print('xpix/ck1 < dt') - - k1 = dt / Cm1 - x1 = (1. - Dm1) / 2. - - # Calc the corrected mass-conservative expression for the reach segment storage at time t+dt - # The lateral inflow ql is only explicitly accounted for in the mass balance equation, while it is not in the equation expressing - # the storage as a weighted average of inflow and outflow.The rationale of this approach lies in the fact that the outflow - # of the reach implicitly takes the effect of the lateral inflow into account. - if q11 == 0: - V11 = V00 + (q00 + q01 - q10 - q11) * dt / 2 - else: - V11 = (1-Dm1)*dt/(2*Cm1)*q01 + (1+Dm1)*dt/(2*Cm1)*q11 - # V11 = k1 * (x1 * q01 + (1. - x1) * q11) # MUST be the same as above! - - ### calc integration on the control volume (pixel) - # calc average discharge outflow q1m for MCT channels during routing sub step dt - # Calculate average outflow using water balance for MCT channel grid cell over sub-routing step - q1mm = q0mm + ql + (V00 - V11) / dt - - # cmcheck - # q1m cannot be smaller than eps or it will cause instability - if q1mm < eps: - q1mm = eps - V11 = V00 + (q0mm + ql - q1mm) * dt - - # q11 Outflow at O(t+dt) - # q1m average outflow in time dt - # V11 water volume at t+dt - # Cm1 Courant number at t+dt for state files - # Dm1 Reynolds numbers at t+dt for state files - return q11, V11, q1mm, Cm1, Dm1 - - - - def hoq(self,q,s0,Balv,ANalv,Nalv): - """Water depth h from discharge q. - Given a generic cross-section (rectangular, triangular or trapezoidal) and a steady-state discharge q=Q*, it computes - water depth (y), wet contour (Bx), wet area (Ax) and wave celerity (cel) using Newton-Raphson method. - Reference: - Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. - Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 - - :param q: steady-state river discharge [m3/s] - :param s0: river bed slope (tan B) - :param Balv : width of the riverbed [m] - :param ChanSdXdY : slope dx/dy of riverbed side - :param ANalv : angle of the riverbed side [rad] - :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] - - :return: - y: water depth referred to the bottom of the riverbed [m] - """ - - alpha = 5./3. # exponent (5/3) - eps = 1.e-06 - max_tries = 1000 - - rs0 = np.sqrt(s0) - usalpha = 1. / alpha - - # cotangent(angle of the riverbed side - dXdY) - if ANalv < np.pi/2: - # triangular or trapezoid cross-section - c = self.cotan(ANalv) - else: - # rectangular corss-section - c = 0. - - # sin(angle of the riverbed side - dXdY) - if ANalv < np.pi/2: - # triangular or trapezoid cross-section - s = np.sin(ANalv) - else: - # rectangular cross-section - s = 1. - - # water depth first approximation y0 based on steady state q - if Balv == 0: - # triangular cross-section - y = (Nalv * q / rs0)**(3. / 8.) * (2 / s)**.25 / c**(5. / 8.) - else: - # rectangular cross-section and first approx for trapezoidal cross-section - y = (Nalv * q / (rs0 * Balv))**usalpha - - if (Balv != 0) and (ANalv < np.pi/2): - # trapezoid cross-section - y = (Nalv * q / rs0) ** usalpha * (Balv + 2. * y / s) ** .4 / (Balv + c * y) - - for tries in range(1,max_tries): - # calc Q(y) for the different tries of y - q0,Ax,Bx,Px,cel = self.qoh(y,s0,Balv,ANalv,Nalv) - # Ax: wet area[m2] - # Bx: cross-section width at water surface[m] - # Px: cross-section wet contour [m] - # cel: wave celerity[m/s] - - # this is the function we want to find the 0 for f(y)=Q(y)-Q* - fy = q0 - q - # calc first derivative of f(y) f'(y)=Bx(y)*cel(y) - dfy = Bx * cel - # calc update for water depth y - dy = fy / dfy - # update yt+1=yt-f'(yt)/f(yt) - y = y - dy - # stop loop if correction becomes too small - if np.abs(dy) < eps: break - - return y - - - def qoh(self,y,s0,Balv,ANalv,Nalv): - """ Discharge q from water depth h. - Given a generic river cross-section (rectangular, triangular and trapezoidal) - and a water depth (y [m]) referred to the bottom of the riverbed, it uses Manning’s formula to calculate: - q: steady-state discharge river discharge [m3/s] - a: wet area [m2] - b: cross-section width at water surface [m] - p: cross-section wet contour [m] - cel: wave celerity [m/s] - Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 - - :param y: river water depth [m] - :param s0: river bed slope (tan B) - :param Balv : width of the riverbed [m] - :param ANalv : angle of the riverbed side [rad] - :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] - - :return: - q,a,b,p,cel - """ - - alpha = 5./3. # exponent (5/3) - - rs0 = np.sqrt(s0) - alpham = alpha - 1. - - # np.where(ANalv < np.pi/2, triang. or trapeiz., rectangular) - # cotangent(angle of the riverbed side - dXdY) - c = np.where(ANalv < np.pi/2, - # triangular or trapezoid cross-section - self.cotan(ANalv), - # rectangular cross-section - 0.) - # sin(angle of the riverbed side - dXdY) - s = np.where(ANalv < np.pi/2, - # triangular or trapezoid cross-section - np.sin(ANalv), - # rectangular corss-section - 1.) - - a = (Balv + y * c) * y # wet area [m2] - b = Balv + 2. * y * c # cross-section width at water surface [m] - p = Balv + 2. * y / s # cross-section wet contour [m] - q = rs0 / Nalv * a**alpha / p**alpham # steady-state discharge [m3/s] - cel = (q / 3.) * (5. / a - 4. / (p * b * s)) # wave celerity [m/s] - - return q,a,b,p,cel - - - def hoV(self,V,xpix,Balv,ANalv): - """Water depth h from volume V. - Given a generic river cross-section (rectangular, triangular and trapezoidal) and a river channel volume V, - it calculates the water depth referred to the bottom of the riverbed [m] (y). - Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 - - :param V : volume of water in channel riverbed - :param xpix : dimension along the flow direction [m] - :param Balv : width of the riverbed [m] - :param ANalv : angle of the riverbed side [rad] - - :return: - y : channel water depth [m] - """ - - - eps = 1e-6 - - c = np.where(ANalv < np.pi/2, # angle of the riverbed side dXdY [rad] - self.cotan(ANalv), # triangular or trapezoidal cross-section - 0.) # rectangular cross-section - - a = V / xpix # wet area [m2] - - # np.where(c < 1.d-6, rectangular, triangular or trapezoidal) - y = np.where(np.abs(c) < eps, - a/Balv, # rectangular cross-section - (-Balv + np.sqrt(Balv**2 + 4 * a * c)) / (2 * c)) # triangular or trapezoidal cross-section - - return y - - - def qoV(self,V,xpix,s0,Balv,ANalv,Nalv): - """ Discharge q from river channel volume V. - Given a generic river cross-section (rectangular, triangular and trapezoidal) - and a water volume (V [m3]), it uses Manning’s formula to calculate the corresponding discharge (q [m3/s]). - Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030 - - :param V : volume of water in channel riverbed - :param xpix : dimension along the flow direction [m] - :param s0: river bed slope (tan B) - :param Balv : width of the riverbed [m] - :param ANalv : angle of the riverbed side [rad] - :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3] - - :return: - y : channel water depth [m] - """ - y = self.hoV(V,xpix,Balv,ANalv) - q, a, b, p, cel = self.qoh(y,s0,Balv,ANalv,Nalv) - return q - - - def cotan(self,x): - """There is no cotangent function in numpy""" - return np.cos(x) / np.sin(x) - - - def get_mct_pix(self,var): + mctpix = self.mct_river_router.pixels_ordered[index] + kinpix = self.river_router.pixels[mapping_mct[mctpix]] + upstream_lookup = self.river_router.upstream_lookup + + # get upstream inflow values from current and previous steps + q00 = 0.0 + q0m = 0.0 + q01 = 0.0 + upstream_pixels = upstream_lookup[kinpix] + for ups_ix in range(upstream_pixels): + ups_pix = upstream_pixels[ups_ix] + q00 += ChanQ_0[ups_pix] + q0m += ChanQAvgDt_0[ups_pix] + q01 += ChanQ_1[ups_pix] + + # get outflow from previous step + q10 = ChanQ_0[kinpix] + V00 = ChanM3_0[kinpix] + ql = SideflowChanMCT[kinpix] + Cm0 = self.var.PrevCm0[kinpix] + Dm0 = self.var.PrevDm0[kinpix] + + # static data + xpix = self.var.ChanLength[kinpix] + s0 = self.var.ChanGrad[kinpix] + Balv = self.var.ChanBottomWidth[kinpix] + Nalv = self.var.ChanManMCT[kinpix] + ANalv = np.arctan(1 / self.var.ChanSdXdY[kinpix]) + + # calling MCT function for single cell idpx + q11, q1m, V11, Cm1, Dm1 = MCTRouting_single( + V00, q10, q01, q00, ql, q0m, Cm0, Dm0, + dt, xpix, s0, Balv, ANalv, Nalv) + + self.var.ChanQ[kinpix] = q11 + self.var.ChanQAvgDt[kinpix] = q1m + self.var.ChanM3[kinpix] = V11 + self.var.PrevCm0[kinpix] = Cm1 + self.var.PrevDm0[kinpix] = Dm1 + + def compress_mct(self, var): """Compress to mct array. For any array (var) with all catchment pixels, it finds the MCT pixels (x) and reduces the dimension of the array (y). @@ -1646,39 +1161,12 @@ def get_mct_pix(self,var): :return: y: same as input array (var) but with only MCT pixels """ - x = np.ma.masked_where(self.var.IsChannelKinematic,var) + x = np.ma.masked_where(self.var.IsChannelKinematic, var) # mask kinematic pixels y = x.compressed() return y - def put_mct_pix(self,var): - """Explode mct array. - For any array (var) with MCT pixels only, it explodes the dimension to all catchment pixels and puts - values from array var in the corresponding MCT pixels. - Pixels not belonging to MCT get value equal to zero. - Uses self.var.IsChannelKinematic to define MCT pixels. - :param var: array (dimension is the number of MCT cells) - :return: - y: same as input array (var) but only all pixels - """ - zeros_array = np.zeros(self.var.IsChannelKinematic.shape) - x = np.ma.masked_where(self.var.IsChannelKinematic, zeros_array) - x[~x.mask] = var - # explode results on the MCT pixels mask (dim=all) - y = x.data - # update results in array (dim=all) - return y - - - def rad_from_dxdy(self,dxdy): - """Calculate radians""" - rad = np.arctan(1 / dxdy) - angle = np.rad2deg(rad) - return rad - - - from .kinematic_wave_parallel import rebuildFlowMatrix, decodeFlowMatrix, streamLookups, topoDistFromSea import pandas as pd class mctWave: