trying to calculate the average discharge for kinematic routing

Author: Cinzia Mazzetti

src/lisflood/hydrological_modules/kinematic_wave_parallel.py

@@ -122,16 +122,21 @@ def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space
         self.kinematic_wave_warning_printed=False
         self.flagnancheck=flagnancheck
         # Parameters for the solution of the discretised Kinematic wave continuity equation
+        self.alpha_channel = alpha_channel
         self.space_delta = space_delta
+        self.inv_time_delta = 1 / time_delta
         self.beta = beta
         self.inv_beta = 1 / beta
         self.b_minus_1 = beta - 1
         self.a_dx_div_dt_channel = alpha_channel * space_delta / time_delta
         self.b_a_dx_div_dt_channel = beta * self.a_dx_div_dt_channel
+        self.a_dx = alpha_channel * space_delta
         # If split-routing (floodplains)
         if alpha_floodplains is not None:
+            self.alpha_floodplains = alpha_floodplains
             self.a_dx_div_dt_floodplains = alpha_floodplains * space_delta / time_delta
             self.b_a_dx_div_dt_floodplains = beta * self.a_dx_div_dt_floodplains
+            self.a_dx_floodplains = alpha_floodplains * space_delta
         # Process flow direction matrix: downstream and upstream lookups, and routing orders
         flow_dir = decodeFlowMatrix(rebuildFlowMatrix(compressed_encoded_ldd, land_mask))
         self.downstream_lookup, self.upstream_lookup = streamLookups(flow_dir, land_mask)
@@ -163,26 +168,28 @@ def _setRoutingOrders(self):
         self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int) 
         self.pixels_ordered = self.pixels_ordered.values.astype(int) 
 
-    def kinematicWaveRouting(self, discharge, specific_lateral_inflow, section="main_channel"):
+    def kinematicWaveRouting(self, discharge_avg, discharge, specific_lateral_inflow, section="main_channel"):
         """Kinematic wave routing: wrapper around kinematic_wave_parallel_tools.kinematicWave"""
         # Lateral inflow (m3 s-1)
         lateral_inflow = nx.evaluate("q * dx", local_dict={"q": specific_lateral_inflow, "dx": self.space_delta})
         # Choose between main channel and floodplain routing
         if section == "main_channel":
             a_dx_div_dt = self.a_dx_div_dt_channel
             b_a_dx_div_dt = self.b_a_dx_div_dt_channel
+            a_dx = self.a_dx
         elif section == "floodplains":
             a_dx_div_dt = self.a_dx_div_dt_floodplains
             b_a_dx_div_dt = self.b_a_dx_div_dt_floodplains
+            a_dx = self.a_dx_floodplains
         else:
             raise Exception("The section parameter must be either 'main_channel' or 'floodplain'!")
         # Constant term in f(x) evaluation for Newton-Raphson method
         local_dict = {"a_dx_div_dt": a_dx_div_dt, "Qold": discharge, "b": self.beta, "lateral_inflow": lateral_inflow}
         constant = nx.evaluate("a_dx_div_dt * Qold ** b + lateral_inflow", local_dict=local_dict)
         # Solve the Kinematic wave equation
-        kwpt.kinematicRouting(discharge, lateral_inflow, constant, self.upstream_lookup,\
+        kwpt.kinematicRouting(discharge_avg, discharge, lateral_inflow, constant, self.upstream_lookup,\
                               self.num_upstream_pixels, self.pixels_ordered, self.order_start_stop,\
-                              self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt)
+                              self.inv_time_delta,self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt, a_dx)
         if self.flagnancheck:
             if self.kinematic_wave_warning_printed==False:
                 if np.all(np.isfinite(discharge))==False:

src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py

@@ -32,9 +32,9 @@
 # ROUTING FUNCTIONS
 # -------------------------------------------------------------------------------------------------
 @njit(parallel=True, fastmath=False, cache=True)
-def kinematicRouting(discharge, lateral_inflow, constant, upstream_lookup,\
-                     num_upstream_pixels, ordered_pixels, start_stop, beta, inv_beta,\
-                     b_minus_1, a_dx_div_dt, b_a_dx_div_dt):
+def kinematicRouting(discharge_avg, discharge, lateral_inflow, constant, upstream_lookup,\
+                     num_upstream_pixels, ordered_pixels, start_stop, inv_time_delta, beta, inv_beta,\
+                     b_minus_1, a_dx_div_dt, b_a_dx_div_dt, a_dx):
     """
     This function performs the kinematic wave routing algorithm to simulate the movement of water through a network of interconnected channels.
     :param discharge:
@@ -49,6 +49,7 @@ def kinematicRouting(discharge, lateral_inflow, constant, upstream_lookup,\
     :param b_minus_1:
     :param a_dx_div_dt:
     :param b_a_dx_div_dt:
+    :param a_dx: ChannelAlpha * ChanLength
     :return:
     """
     num_orders = start_stop.shape[0]
@@ -59,34 +60,43 @@ def kinematicRouting(discharge, lateral_inflow, constant, upstream_lookup,\
         # Iterate through each pixel in the current order in parallel
         for index in prange(first, last):
             # Solve the kinematic wave for the current pixel
-            solve1Pixel(ordered_pixels[index], discharge, lateral_inflow, constant, upstream_lookup,\
-                        num_upstream_pixels, a_dx_div_dt, b_a_dx_div_dt, beta, inv_beta, b_minus_1)
+            solve1Pixel(ordered_pixels[index], discharge_avg, discharge, lateral_inflow, constant, upstream_lookup,\
+                        num_upstream_pixels, a_dx_div_dt, b_a_dx_div_dt, inv_time_delta, beta, inv_beta, b_minus_1, a_dx)
 
 @njit(nogil=True, fastmath=False, cache=True)
-def solve1Pixel(pix, discharge, lateral_inflow, constant,\
+def solve1Pixel(pix, discharge_avg, discharge, lateral_inflow, constant,\
                       upstream_lookup, num_upstream_pixels, a_dx_div_dt,\
-                      b_a_dx_div_dt, beta, inv_beta, b_minus_1):
+                      b_a_dx_div_dt, inv_time_delta, beta, inv_beta, b_minus_1, a_dx):
     """
     Te Chow et al. 1988 - Applied Hydrology - Chapter 9.6
     :param pix:
-    :param discharge:
+    :param discharge_avg: average outflow discharge
+    :param discharge: instantaneous outflow discharge
     :param lateral_inflow:
     :param constant:
     :param upstream_lookup:
     :param num_upstream_pixels:
     :param a_dx_div_dt:
     :param b_a_dx_div_dt:
+    :param inv_time_delta: 1/DtRouting
     :param beta:
     :param inv_beta:
     :param b_minus_1:
+    :param a_dx: ChannelAlpha * ChanLength
     :return:
     """
     count = 0
     previous_estimate = -1.0
     upstream_inflow = 0.0
+    upstream_inflow_avg = 0.0
+
+    # volume of water in channel at beginning of computation step
+    channel_volume_start = a_dx * discharge[pix]**beta
+
     # Inflow from upstream pixels
     for ups_ix in range(num_upstream_pixels[pix]):
         upstream_inflow += discharge[upstream_lookup[pix,ups_ix]]
+        upstream_inflow_avg += discharge_avg[upstream_lookup[pix, ups_ix]]
     const_plus_ups_infl = upstream_inflow + constant[pix] # upstream_inflow + alpha*dx/dt*Qold**beta + dx*specific_lateral_inflow
     # If old discharge, upstream inflow and lateral inflow are below accuracy: set discharge to 0 and exit
     if const_plus_ups_infl <= NEWTON_TOL:
@@ -111,6 +121,17 @@ def solve1Pixel(pix, discharge, lateral_inflow, constant,\
     # If iterations converge to NEWTON_TOL, set value to 0
     if discharge[pix] == NEWTON_TOL:
         discharge[pix] = 0
+
+    # cmcheck
+    # avoid negative discharge
+    if discharge[pix] < 0:
+        discharge[pix] = 0
+    # volume of water in channel at end of computation step
+    channel_volume_end = a_dx * discharge[pix]**beta
+    # mass water balance to calc average outflow
+    discharge_avg[pix] = upstream_inflow_avg + lateral_inflow + (channel_volume_start - channel_volume_end) * inv_time_delta
+
+
     # to simulate inf or nan: discharge[pix] = 1.0/0.0
     # with gil:
     #    got_valid_value = np.isfinite(discharge[pix])