Feature/mct restore calib

VERSION

@@ -1 +1 @@
-4.3.1.505
+4.3.1.506

src/lisflood/Lisflood_initial.py

@@ -138,6 +138,7 @@ def __init__(self):
         self.surface_routing_module = surface_routing(self)
         self.reservoir_module = Reservoir(self) # get_reservoir(option['reservoirHanazaki'])
         self.lakes_module = lakes(self)
+
         self.polder_module = polder(self)
         self.waterabstraction_module = waterabstraction(self)
         self.indicatorcalc_module = indicatorcalc(self)
@@ -189,10 +190,12 @@ def __init__(self):
 
         self.snow_module.initial()
         self.frost_module.initial()
+
         self.leafarea_module.initial()
-        self.soilloop_module.initial()
 
+        self.soilloop_module.initial()
         self.soil_module.initial()
+
         self.routing_module.initial()
 
         self.groundwater_module.initial()
@@ -201,6 +204,8 @@ def __init__(self):
         self.inflow_module.initial()
         self.surface_routing_module.initial()
 
+        # At this point LddChan and LddKinematic do not have any structure reservoirs/lakes MCT confluence
+
         self.reservoir_module.initial()
         self.lakes_module.initial()
         self.polder_module.initial()
@@ -209,18 +214,30 @@ def __init__(self):
 
         self.structures_module.initial()
         # Structures such as reservoirs and lakes are modelled by interrupting the channel flow paths
+        # At this point LddKinematic and LddChan have pits upstream of (structures) reservoirs and lakes
 
         # ----------------------------------------------------------------------
         # ----------------------------------------------------------------------
 
         self.routing_module.initialSecond()
         # CHANNEL INITIAL SPLIT UP IN SECOND CHANNEL
+
         self.surface_routing_module.initialSecond()
 
+        # MCT confluence must be in the LddKinematic when I get here
+        self.routing_module.initialKinematicWave()
+
         if option.get('MCTRouting'):
             self.routing_module.initialMCT()
             # initialising Muskingum-Cunge-Todini routing for channel
 
+        # self.routing_module.initialKinematicWave()
+        # this cannot be here because I need river_router in the MCT initialization
+
+
+
+
+
         self.evapowater_module.initial()
         self.riceirrigation_module.initial()
         self.waterabstraction_module.initial()
@@ -397,4 +414,4 @@ def defsoil(self, name_1, name_2=None, name_3=None, coords=None):
     def deffraction(self, variable):
         """Weighted sum over the soil fractions of each pixel"""
         ax_veg = variable.dims.index("vegetation")
-        return _vegSum(ax_veg, variable.values, self.SoilFraction.values)
+        return _vegSum(ax_veg, variable.values, self.SoilFraction.values)

src/lisflood/global_modules/default_options.py

@@ -1288,6 +1288,7 @@
                    'repwateruseGauges': False,
                    'repwateruseSites': False,
                    'riceIrrigation': False,
+                   'simulateCalibrationPoints': False,
                    'simulateLakes': False,
                    'simulatePF': False,
                    'simulatePolders': False,

src/lisflood/hydrological_modules/mct.py

@@ -20,11 +20,13 @@ def __init__(
             dt,                         # computation time step for routing [s]
             river_router,               # class
             mapping_mct,                # MCT pixels mapping
+            CalibPointsIds,             # calibration points pixels ids
         ):
 
         # 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)
+        # Inside streamLookups each land_mask pixel is assigned a unique id from 0 to numpix-1, numbered by row
         self.num_upstream_pixels = (self.upstream_lookup != -1).sum(1).astype(int) # astype for cython import in windows (to avoid 'long long' buffer dtype mismatch)
         # Routing order: decompose domain into batches; within each batch, pixels can be routed in parallel
         self._setMCTRoutingOrders()
@@ -37,7 +39,7 @@ def __init__(
         self.dt = dt
         self.river_router = river_router
         self.mapping_mct = mapping_mct
-
+        self.CalibPointsIds = CalibPointsIds
 
     def _setMCTRoutingOrders(self):
         """Compute the MCT wave routing order. Pixels are grouped in sets with the same order.
@@ -103,6 +105,7 @@ def routing(
             PrevCm0,            # Courant number in input: at time t; in output: at time t+dt
             PrevDm0,            # Reynolds number in input: at time t; in output: at time t+dt
             ChanM3,             # V11 as output
+            self.CalibPointsIds,# inflow points used by the calibration suite
         )
 
 
@@ -130,6 +133,7 @@ def mct_routing(
     PrevCm0,        # Courant number in input: at time t; in output: at time t+dt
     PrevDm0,        # Reynolds number in input: at time t; in output: at time t+dt
     ChanM3,         # V11 as output
+    CalibPointsIds, # inflow points used by the calibration suite
 ):
     """This function implements Muskingum-Cunge-Todini routing method
     MCT routing is calculated on MCT pixels only but gets inflow from both Kinematic/Split and MCT upstream pixels.
@@ -167,11 +171,22 @@ def mct_routing(
             q00 = 0.0
             q0m = 0.0
             q01 = 0.0
+            ql = SideflowChanMCT[kinpix]   # ← ql defined BEFORE the loop
             for ups_ix in range(num_upstream_pixels[kinpix]):
                 ups_pix = upstream_pixels[ups_ix]   # upstream pixel id
-                q00 += ChanQ_0[ups_pix]     # Inflow (x) to the pixel at previous step t (instant)
-                q0m += ChanQAvgDt[ups_pix]  # Average inflow (x) to the pixel at previous step t (average)
-                q01 += ChanQ[ups_pix]       # Inflow (x) at current step t+dt (instant)
+                #####################################################################################################
+                # This is necessary for EFAS6/GloFAs5 calibration
+                if np.any(CalibPointsIds == ups_pix):
+                    # this upstream pixel is a calibration point - add to sideflow
+                    ql += ChanQAvgDt[ups_pix]   # avoid += 
+                    # Sideflow during step dt including contribution from calibration pixel
+                else:
+                    # not a calibration point - go as usual
+                    q00 += ChanQ_0[ups_pix]  # Inflow (x) to the pixel at previous step t (instant)
+                    q0m += ChanQAvgDt[ups_pix]  # Average inflow (x) to the pixel at previous step t (average)
+                    q01 += ChanQ[ups_pix]  # Inflow (x) at current step t+dt (instant)
+
+                #####################################################################################################
 
             # get outflow from the pixel at previous step t
             q10 = ChanQ_0[kinpix]   # Outflow (x+dx) from the pixel at previous step t (instant)
@@ -181,7 +196,7 @@ def mct_routing(
             Cm0 = PrevCm0[kinpix]   # Courant number at the end of previous step t
             Dm0 = PrevDm0[kinpix]   # Reynolds number at the end of previous step t
 
-            ql = SideflowChanMCT[kinpix]    # Sideflow during step dt
+            # ql = SideflowChanMCT[kinpix]    # Sideflow during step dt
 
             # static data
             xpix = ChanLength[kinpix]                   # Channel length
@@ -239,7 +254,7 @@ def MCTRouting_single(
 
     # 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)
+    q11 = q10 + (q01 - q00) # + ql
 
     # check for negative and zero discharge values
     # zero outflow is not allowed
@@ -250,6 +265,11 @@ def MCTRouting_single(
     # qm0 = (I(t)+O(t))/2
     # qm0 = (q00 + q10) / 2.
 
+    # ql correction for hydraulic state computation
+    # Only positive ql contributes (ignore abstraction), 
+    # /2 represents uniform distribution along reach
+    # ql_correction = max(ql, 0.0) / 2.0
+
     # Calc O(t+dt)=q11 at time t+dt using MCT equations
     for i in range(2):  # repeat 2 times for accuracy
 
@@ -273,7 +293,7 @@ def MCTRouting_single(
 
         # Calc reference discharge time t+dt
         # Q(t+dt)=(I(t+dt)+O'(t+dt))/2
-        qm1 = (q01 + q11) / 2.0
+        qm1 = (q01 + q11) / 2.0 # + ql_correction
         # cm
         if qm1 <= eps :  # cmcheck ==0     #tpk
             qm1 = eps                   #tpk
@@ -291,10 +311,24 @@ def MCTRouting_single(
         Cm1 = ck1 * dt / xpix / Beta1
 
         # Calc MCT parameters
+        # Guard Diffusivity
+        # Dm1 = min(max(Dm1, 0.0), 1.0)
         den = 1 + Cm1 + Dm1
+
         c1 = (-1 + Cm1 + Dm1) / den
-        c2 = (1 + Cm0 - Dm0) / den * (Cm1 / Cm0)
-        c3 = (1 - Cm0 + Dm0) / den * (Cm1 / Cm0)
+
+        # Guard against Cm0 near zero (first timestep or after dry conditions)
+        if Cm0 > eps:
+            cm_ratio = Cm1 / Cm0
+            if cm_ratio > 3.0:
+                cm_ratio = 3.0
+            elif cm_ratio < 0.1:
+                cm_ratio = 0.1
+        else:
+            cm_ratio = 1.0
+
+        c2 = (1 + Cm0 - Dm0) / den * cm_ratio
+        c3 = (1 - Cm0 + Dm0) / den * cm_ratio
         c4 = (2 * Cm1) / den
 
         # cmcheck
@@ -304,9 +338,23 @@ def MCTRouting_single(
         # Mass balance equation that takes into consideration the lateral flow
         q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql
 
+        # Diagnostic - now only fires when the CLAMPED ratio is at the ceiling
+        # or other suspicious conditions
+        # if cm_ratio >= 3.0 - 1e-9:  # ratio was clamped
+        #     print("CLAMPED Cm0=", Cm0, "Cm1=", Cm1, "raw_ratio=", Cm1/max(Cm0, 0.001), "q01=", q01)
+
         if q11 < 0:  # cmcheck <=0  #tpk
             q11 = 0                 #tpk
 
+        # NEW: mass-balance upper bound
+        # if ql > 0:
+        #     ql_abs = ql
+        # else:
+        #     ql_abs = -ql
+        max_q11 = q01 + q00 + abs(ql) + V00 / dt
+        if q11 > max_q11:
+            q11 = max_q11
+
         #### end of for loop
 
     # # cmcheck
@@ -326,6 +374,10 @@ def MCTRouting_single(
     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!
+    # transition = min(1.0, q11 / (q11 + eps))  # Sigmoid-like transition
+    # V11_formula1 = V00 + (q00 + q01 - q10 - q11) * dt / 2
+    # V11_formula2 = (1 - Dm1) * dt / (2 * Cm1) * q01 + (1 + Dm1) * dt / (2 * Cm1) * q11
+    # V11 = transition * V11_formula2 + (1 - transition) * V11_formula1
 
     if V11 < 0 :    #tpk
         V11 = 0     #tpk
@@ -334,6 +386,14 @@ def MCTRouting_single(
     # 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
+    # Ensure q1mm is consistent with q11 (instantaneous outflow)
+    # q1mm should be within reasonable bounds of q11
+    # if q11 > 0:
+    #     ratio = q1mm / q11
+    #     if ratio > 10:
+    #         q1mm = q11
+    #     elif ratio < 0.1:
+    #         q1mm = q11
 
     # cmcheck
     # q1m cannot be smaller than eps or it will cause instability