From eb2eb8cc259eaca3080914b5b6ca9e4473b63125 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Mon, 18 May 2026 13:24:16 +0200 Subject: [PATCH 01/12] Own clamp + Cinzias mctconfluence.py --- src/lisflood/hydrological_modules/mct.py | 59 +++++- .../hydrological_modules/mctconfluence.py | 171 ++++++++++++++++++ src/lisflood/hydrological_modules/routing.py | 45 ++++- 3 files changed, 265 insertions(+), 10 deletions(-) create mode 100644 src/lisflood/hydrological_modules/mctconfluence.py diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index 99f0ad22..7490564c 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/src/lisflood/hydrological_modules/mct.py @@ -239,7 +239,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 +250,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 +278,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 +296,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 +323,29 @@ 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 + ql_abs + V00 / dt + # if q11 > max_q11: + # q11 = max_q11 + + # After q11 is computed + # max_physical = q01 + q00 + abs(ql) + V00 / dt + # if q11 > max_physical * 1.1: # exceeds physical bound by >10% + # print("MASS VIOLATION q11=", q11, "max_phys=", max_physical, + # "q01=", q01, "q00=", q00, "ql=", ql, "V00=", V00) + #### end of for loop # # cmcheck @@ -326,6 +365,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 +377,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 = 10 * q11 + # elif ratio < 0.1: + # q1mm = 0.1 * q11 # cmcheck # q1m cannot be smaller than eps or it will cause instability diff --git a/src/lisflood/hydrological_modules/mctconfluence.py b/src/lisflood/hydrological_modules/mctconfluence.py new file mode 100644 index 00000000..fa670825 --- /dev/null +++ b/src/lisflood/hydrological_modules/mctconfluence.py @@ -0,0 +1,171 @@ +""" + +Copyright 2019 European Union + +Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission +subsequent versions of the EUPL (the "Licence"); + +You may not use this work except in compliance with the Licence. +You may obtain a copy of the Licence at: + +https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt + +Unless required by applicable law or agreed to in writing, +software distributed under the Licence is distributed on an "AS IS" basis, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the Licence for the specific language governing permissions and limitations under the Licence. + +""" + +from __future__ import print_function, absolute_import +from pcraster import Scalar, numpy2pcr, pcr2numpy,downstream, boolean,lddrepair,ifthenelse,lddmask +from nine import range +import numpy as np +from ..global_modules.settings import LisSettings, MaskInfo +from ..global_modules.add1 import loadmap, compressArray, decompress, makenumpy +from ..global_modules.errors import LisfloodWarning +from . import HydroModule + + +class mctconfluence(HydroModule): + """ + Adds contribution from Kinematic cells to MCT cells as a lateral flow when MCT routing is enabled. + This is for MCT cells that have upstream contributions from both kinematic and MCT cells. + + This module handles the initialization and dynamic simulation of Kinematic to MCT cells confluence. + It injects side discharge to the downstream grid cell as lateral flow. + + Attributes: + ----------- + var (object): An object containing all the variables used within the mctconfluence module. + + Methods: + -------- + initial(): Sets up the initial conditions and parameters for the simulation, + including confluence locations. + dynamic_inloop(NoRoutingExecuted: int): Performs dynamic calculations within the routing + loop to simulate the kinematic to MCT confluence. + """ + + module_name = 'MCTConfluence' + + def __init__(self, mctconfluence_variable): + self.var = mctconfluence_variable + + def __init__(self, mctconfluence_variable): + """ + Initializes the MCT confluence module with a given variable object. + + Parameters: + ----------- + mctconfluence_variable: object + An object containing the variables needed for the MCT confluence simulation. + """ + + self.var = mctconfluence_variable + + def initial(self): + """ + Initiates the MCT confluence module by loading the necessary data and maps. + """ + + settings = LisSettings.instance() + option = settings.options + binding = settings.binding + maskinfo = MaskInfo.instance() + if option['MCTRouting']: + + inArPcr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32")) # pcr + # Assign a number to each non-missing pixel as cell id, starting from 0 + inAr = compressArray(inArPcr) #np + + down = (compressArray(downstream(self.var.LddStructuresChan, inArPcr))).astype("int32") # np + # assign to each pixel the cell id of the pixel it is contributing to + # LddStructuresChan do not contain structures, LddStructuresKinematic is same as LddStructuresChan + + LddKinematic = lddmask(self.var.LddStructuresChan, self.var.IsChannelKinematicPcr) + # Mask of LddKinematic with only kinematic cells and no structures + # Sinks are added at the last Kin pixel before confluence with MCT pixels + + maskKinematic = (compressArray(LddKinematic) == 5) #& (self.var.IsUpsOfStructureKinematicC != 1) + + # find location of KIN pixels (only) in LddKin that are upstream of the confluence with an MCT pixel + maskKinematic[compressArray(self.var.AtLastPoint) == 1] = False + # remove sinks that are outlets + # this does NOT include sinks upstream of structures (lakes, reservoirs) and outlets + maskKinematicPcr = boolean(decompress(maskKinematic)) + + self.var.LddChan = ifthenelse(maskKinematicPcr, 5, self.var.LddChan) + self.var.LddKinematic = ifthenelse(maskKinematicPcr, 5, self.var.LddKinematic) + # Adding sinks to Ldd at the last KIN pixels upstream of the confluence with an MCT pixel to LddChan and LddKinematic + # This similar to what is done in structures + # At this point LddChan and LddKinematic have sinks upstteam of structures and of a Kin-MCT confluence and outlets + + self.var.KinematicUpsOfMCTConfluence = np.where(maskKinematic, down, 0) + # find last KIN pixels upstream of the confluence with an MCT pixel and assign it the id of the downstream MCT pixel + + mctconfluence = maskinfo.in_zero() + mctconfluence[np.isin(inAr, self.var.KinematicUpsOfMCTConfluence[self.var.KinematicUpsOfMCTConfluence != 0])] = 1 + # for each element in KinematicUpsOfMCTConfluence (they are Kinematic cells), get the value of the downstream cell, + # then find the position ix of that same value in inAr (vector with numbering of all cells) this is the position of the MCT confluence cell, + # read the inAr value and put 1 in the corrisponding position in mctconfluence + # identify location of MCT pixels that receive a contribution from an upstream KIN pixel + + mctconfluence = np.where(self.var.IsStructureChan, 0, mctconfluence) + # remove the MCT cells with a structure (reservoir/lake) on them + + self.var.MCTConfluenceSitesC = mctconfluence + self.var.MCTConfluenceSitesCC = np.compress(mctconfluence > 0, mctconfluence) + self.var.MCTConfluenceIndex = np.nonzero(mctconfluence)[0] + pass + + + def dynamic_init(self): + """ Initialization of the dynamic part of the MCT confluence module + init mct confluence before sub step routing + """ + settings = LisSettings.instance() + option = settings.options + if option['MCTRouting']: + self.var.QInConfM3Old = np.where(self.var.MCTConfluenceSitesC > 0, self.var.ChanQAvgDt * self.var.DtSec, 0) + + + def dynamic_inloop(self, NoRoutingExecuted: int): + """ + Performs the dynamic simulation of MCT confluence within the routing loop. This method + injects upstream discharge to the downstream grid cell as lateral inflow. + + Parameters: + ----------- + NoRoutingExecuted: integer + The number of routing sub-steps that have been executed. This parameter is used to manage + the accumulation of inflow and outflow over the routing steps. + """ + + settings = LisSettings.instance() + option = settings.options + maskinfo = MaskInfo.instance() + + # self.var.QConfADDEDM3 = maskinfo.in_zero() + + if option['MCTRouting'] and not option['InitLisflood']: + + InvDtSecDay = 1 / float(86400) + # InvDtSecDay=self.var.InvDtSec + + lateralflow = np.bincount(self.var.KinematicUpsOfMCTConfluence, weights=self.var.ChanQAvgDt)[self.var.MCTConfluenceIndex] #same as Qin + # contribution to the MCT pixel from upstream Kinematic pixels + + self.var.QInConfM3 = maskinfo.in_zero() + np.put(self.var.QInConfM3, self.var.MCTConfluenceIndex, lateralflow * self.var.DtSec) + self.var.QDeltaConfM3 = (self.var.QInConfM3 - self.var.QInConfM3Old) * self.var.InvNoRoutSteps + # difference between old and new lateral flow per sub step + # in order to calculate the amount of lateral flow in the routing loop + + self.var.QConfM3Dt = (self.var.QInConfM3Old + (NoRoutingExecuted + 1) * self.var.QDeltaConfM3) * self.var.InvNoRoutSteps + # output to the MCT confluence cell + + self.var.QInConfM3Old = self.var.QInConfM3.copy() + # save the lateral flow for next step + + # self.var.QConfADDEDM3 += self.var.QConfM3Dt \ No newline at end of file diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index feb2fc7d..ee610db9 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -459,7 +459,7 @@ def initialSecond(self): # Over bankful discharge starts at QLimit # lower discharge limit for second line of routing # set to multiple of average discharge (map from prerun) - # QSplitMult =2 is around 90 to 95% of Q + # QSplitMult = 2 is around 90 to 95% of Q self.var.QLimit = loadmap_base('AvgDis') * loadmap('QSplitMult') # Using loadmap_base function as we don't want to cache avgdis in the calibration @@ -807,7 +807,32 @@ def dynamic(self, NoRoutingExecuted): self.var.PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt self.var.ChanM3 # Channel storage volume. In input: at time t V00; in output: at time t+dt V11 ) + problematic = np.where(self.var.PrevDm0[self.var.IsChannelMCT] > 2)[0] + # Case 1: ChanQAvgDt too small relative to ChanQ (flood arrival) + mct_avg_too_small = ( + self.var.IsChannelMCT & + (self.var.ChanQ > 100.) & + (self.var.ChanQAvgDt < 0.1 * self.var.ChanQ) + ) + self.var.ChanQAvgDt = np.where(mct_avg_too_small, self.var.ChanQ, self.var.ChanQAvgDt) + + # Case 2: ChanQAvgDt too large relative to ChanQ (flood recession / cold start) + mct_avg_too_large = ( + self.var.IsChannelMCT & + (self.var.ChanQAvgDt > 100.) & + (self.var.ChanQ < 0.1 * self.var.ChanQAvgDt) + ) + self.var.ChanQAvgDt = np.where(mct_avg_too_large, self.var.ChanQ, self.var.ChanQAvgDt) + # if 95 <= self.var.currentStep <= 110: + # mct_pixels = self.var.IsChannelMCT + # cm0 = self.var.PrevCm0[mct_pixels] # this is now Cm1 after routing + # # We can't get cm_ratio directly but we can see Cm distribution + # print(f"Step {self.var.currentStep}, substep {NoRoutingExecuted}:") + # print(f" Cm: min={cm0.min():.2f} mean={cm0.mean():.2f} max={cm0.max():.2f}") + # print(f" Dm: min={self.var.PrevDm0[mct_pixels].min():.2f} " + # f"mean={self.var.PrevDm0[mct_pixels].mean():.2f} " + # f"max={self.var.PrevDm0[mct_pixels].max():.2f}") ##################################################################3 if flags['debug']: # checking Courant number for potential instability in MCT @@ -824,11 +849,19 @@ def dynamic(self, NoRoutingExecuted): bad = too_large | too_small - if np.any(bad): - warnings.warn(LisfloodWarning("WARNING! At least one ChanQ is >> or << ChanQAvgDt. Consider increasing DtRouting step or using kinematic routing")) - # # list 'bad' cells - # bad_indices = np.where(dismask)[0][bad] - # print("Bad indices:", bad_indices) + # if np.any(bad): + # warnings.warn(LisfloodWarning("WARNING! At least one ChanQ is >> or << ChanQAvgDt. Consider increasing DtRouting step or using kinematic routing")) + # # # list 'bad' cells + # # bad_indices = np.where(dismask)[0][bad] + # # print("Bad indices:", bad_indices) + # bad_indices = np.where(dismask)[0][bad] + # print("Step:", NoRoutingExecuted) + # print("Bad indices:", bad_indices) + # print("ChanQ at bad:", self.var.ChanQ[bad_indices]) + # print("ChanQAvgDt at bad:", self.var.ChanQAvgDt[bad_indices]) + # print("PrevCm0 at bad:", self.var.PrevCm0[bad_indices]) + # print("PrevDm0 at bad:", self.var.PrevDm0[bad_indices]) + # print("IsChannelMCT at bad:", self.var.IsChannelMCT[bad_indices]) ##################################################################3 else: From 723532d8ba72e2c85b1bfbdf52029793ca81b1d3 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Mon, 18 May 2026 15:40:32 +0200 Subject: [PATCH 02/12] Further Cinzia mct corrections --- src/lisflood/Lisflood_initial.py | 31 +++++- src/lisflood/hydrological_modules/routing.py | 106 ++++++++++++++----- 2 files changed, 108 insertions(+), 29 deletions(-) diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index 346c1dff..084377dd 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/src/lisflood/Lisflood_initial.py @@ -49,6 +49,7 @@ from .hydrological_modules.surface_routing import surface_routing from .hydrological_modules.reservoir import Reservoir from .hydrological_modules.lakes import lakes +from .hydrological_modules.mctconfluence import mctconfluence from .hydrological_modules.polder import polder from .hydrological_modules.waterabstraction import waterabstraction from .hydrological_modules.indicatorcalc import indicatorcalc @@ -138,6 +139,8 @@ 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.mctheadwater_module = mctheadwater(self) + self.mctconfluence_module = mctconfluence(self) self.polder_module = polder(self) self.waterabstraction_module = waterabstraction(self) self.indicatorcalc_module = indicatorcalc(self) @@ -189,10 +192,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 +206,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 +216,38 @@ 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.mctheadwater_module.initial() + # # initialising MCT headwater pixels and adding pixels upstream of MCT headwater to the LDD + self.mctconfluence_module.initial() + # initialising MCT confluence points and adding MCT confluence sinks to the LDD # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- 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.mctheadwater_module.dynamic_init() + self.mctconfluence_module.dynamic_init() + + # 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 +424,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) \ No newline at end of file diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index ee610db9..2adc1074 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -16,7 +16,7 @@ """ from __future__ import print_function, absolute_import -from pcraster import lddmask, accuflux, boolean, downstream, pit, path, lddrepair, ifthenelse, cover, nominal, uniqueid, \ +from pcraster import lddmask, accuflux, boolean, scalar, downstream, pit, path, lddrepair, ifthenelse, cover, nominal, uniqueid, \ catchment, upstream, pcr2numpy import warnings @@ -30,6 +30,7 @@ from .transmission import transmission from .kinematic_wave_parallel import kinematicWave, kwpt from .mct import MCTWave +from .mctconfluence import mctconfluence from ..global_modules.settings import LisSettings, MaskInfo from ..global_modules.errors import LisfloodWarning @@ -61,6 +62,7 @@ def __init__(self, routing_variable): self.polder_module = polder(self.var) self.inflow_module = inflow(self.var) self.transmission_module = transmission(self.var) + self.mctconfluence_module = mctconfluence(self.var) # -------------------------------------------------------------------------- # -------------------------------------------------------------------------- @@ -140,6 +142,7 @@ def initial(self): self.var.LddToChan = lddrepair(ifthenelse(self.var.IsChannelPcr, 5, self.var.Ldd)) #pcr self.var.LddToChanNp=compressArray(self.var.LddToChan) #np # Routing of runoff (incl. groundwater) to the river channel + # LDD for routing runoff (incl. groundwater) to the channel if option['dynamicWave']: pass @@ -174,6 +177,45 @@ def initial(self): self.var.LddKinematic = self.var.LddChan self.var.LddKinematicNp = compressArray(self.var.LddKinematic) # np + + # At this point, LddChan and LddKinematic do not have sinks at reservoirs/lakes or MCT headwater and MCT confluences + # LddMCT does not exist yet + + # ************************************************************ + # ***** MCT DRAINAGE NETWORK GEOMETRY - LDD ***************** + # ************************************************************ + + # This is done here to be able to add MCT headwater pixels to structures + if option['MCTRouting']: + + self.var.IsChannelMCTPcr = boolean(loadmap('ChannelsMCT', pcr=True)) # pcr + # load mask of MCT river grid cells + self.var.IsChannelMCT = np.bool8(compressArray(self.var.IsChannelMCTPcr)) # bool + + # even if MCT is active, it should be deactivated if there is no MCT cell in the domain + if self.var.IsChannelMCT.sum() == 0: + warnings.warn(LisfloodWarning('There are no MCT grid cell. MCT routing is deactivated')) + option['MCTRouting'] = False + # rebuild lists of reported files with MCTRouting = False + settings.build_reportedmaps_dicts() + + if option['MCTRouting'] and not option['InitLisflood']: + + self.var.IsChannelMCTPcr = boolean(decompress(self.var.IsChannelMCT)) # pcr + # Identify channel pixels where Muskingum-Cunge-Todini is used + + self.var.mctmask = np.bool8(pcr2numpy(self.var.IsChannelMCTPcr,0)) + # Create a mask with cells using MCT + + self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr + self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np + # Identify channel pixels where Kinematic wave is used instead of MCT + + + # ************************************************************ + # ***** MCT DRAINAGE NETWORK GEOMETRY - LDD ***************** + # ************************************************************ + self.var.AtLastPoint = boolean(pit(self.var.Ldd)) #pcr # Assign True to each of the grid cells where there are outlet points # Function 'pit' assigns a unique number starting from 1 to pit cells (ldd=5) in the Ldd @@ -248,6 +290,10 @@ def initial(self): TotalCrossSectionAreaHalfBankFull = BankFullPerc * self.var.TotalCrossSectionAreaBankFull # set BankFullPerc to 0.5 for half bankfull + # Channel volume initialization for MCT cells + TotalCrossSectionAreaHalfBankFull = np.where(self.var.IsChannelKinematic, TotalCrossSectionAreaHalfBankFull, 0.01 * self.var.TotalCrossSectionAreaBankFull) + # set initial volume in MCT cells to 1% of bankfull + TotalCrossSectionAreaInitValue = loadmap('TotalCrossSectionAreaInitValue') self.var.TotalCrossSectionArea = np.where(TotalCrossSectionAreaInitValue == -9999, TotalCrossSectionAreaHalfBankFull, TotalCrossSectionAreaInitValue) # Total cross-sectional area [m2]: if initial value in binding equals -9999 the value at half bankfull is used, @@ -486,7 +532,15 @@ def initialSecond(self): # Total (virtual) outflow from river channel when second routing line is active (= using increased Manning coeff 2) self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta) # (Real) outflow from main channel when second line of routing is active (= using riverbed Manning coeff 2) - + + def initialKinematicWave(self): + """ Initialization of the parallel kinematic wave router for Kinematic routing and SplitRouting: + main channel-only routing if self.var.ChannelAlpha2 is None; else split-routing(main channel + floodplains). + Initialization uses LDD for kinematic routing (LddKinematic) + """ + settings = LisSettings.instance() + option = settings.options + flags = settings.flags # ************************************************************ # ***** INITIALISE PARALLEL KINEMATIC WAVE ROUTER ************ @@ -505,7 +559,7 @@ def initialSecond(self): if option['InitLisflood'] and option['repMBTs']: # Calculate initial water storage in rivers (no lakes no reservoirs) # self.var.StorageStepINIT= self.var.ChanM3Kin - self.var.StorageStepINIT = self.var.ChanM3 + self.var.StorageStepINIT = self.var.ChanM3.copy() # Initial water volume in river channels self.var.DischargeM3StructuresIni = maskinfo.in_zero() if option['simulateReservoirs']: @@ -515,7 +569,7 @@ def initialSecond(self): self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments) if not option['InitLisflood'] and option['repMBTs']: - self.var.StorageStepINIT = self.var.ChanM3 + self.var.StorageStepINIT = self.var.ChanM3.copy() # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0) DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0) if not(option['SplitRouting']): @@ -554,28 +608,13 @@ def initialMCT(self): # ***** INITIALISATION FOR MCT ROUTING ************ # ************************************************************ - # even if MCT is active, it should be deactivated if there is no MCT cell in the domain - if option['MCTRouting']: - self.var.IsChannelMCTPcr = boolean(loadmap('ChannelsMCT', pcr=True)) #pcr - self.var.IsChannelMCT = np.bool8(compressArray(self.var.IsChannelMCTPcr)) #bool - if self.var.IsChannelMCT.sum()==0: - warnings.warn(LisfloodWarning('There are no MCT grid cell. MCT routing is deactivated')) - option['MCTRouting'] = False - # rebuild lists of reported files with MCTRouting = False - settings.build_reportedmaps_dicts() if option['MCTRouting'] and not option['InitLisflood']: maskinfo = MaskInfo.instance() - self.var.IsChannelMCTPcr = boolean(decompress(self.var.IsChannelMCT)) # pcr - # Identify channel pixels where Muskingum-Cunge-Todini is used - - self.var.mctmask = np.bool8(pcr2numpy(self.var.IsChannelMCTPcr,0)) - # Create a mask with cells using MCT - - self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr - self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np - # Identify channel pixels where Kinematic wave is used instead of MCT + # self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr + # self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np + # # Identify channel pixels where Kinematic wave is used instead of MCT self.var.LddMCT = lddmask(self.var.LddChan, self.var.IsChannelMCTPcr) #pcr # Ldd for MCT routing @@ -591,7 +630,6 @@ def initialMCT(self): self.var.ChanGrad[MCT_slope_mask] = ChanGradMaxMCT # set max channel slope for MCT pixels - # cmcheck # This could become a calibration parameter if we want to use MCT+SplitRouting self.var.ChanManMCT = (self.var.ChanMan / self.var.CalChanMan) * loadmap('CalChanMan3') # Mannings coefficient for MCT pixels (same as second line of split routing) @@ -603,6 +641,7 @@ def initialMCT(self): self.var.PrevDm0 = np.where(PrevDmMCT == -9999, maskinfo.in_zero(), PrevDmMCT) #np # Reynolds number (Dm) for MCT at previous time step t0 + # ************************************************************ # ***** INITIALISE MUSKINGUM-CUNGE-TODINI WAVE ROUTER ******** # ************************************************************ @@ -702,7 +741,7 @@ def dynamic(self, NoRoutingExecuted): else: self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.ToChanM3RunoffDt.copy()),self.var.Catchments) if option['inflow']: - self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.QInDt),self.var.Catchments) + self.var.AddedTRUN += np.take(np.bincount(self.var.Catchments, weights=self.var.QInDt),self.var.Catchments) if option['openwaterevapo']: self.var.AddedTRUN -= np.take(np.bincount(self.var.Catchments, weights=self.var.EvaAddM3Dt.copy()),self.var.Catchments) if option['wateruse']: @@ -782,9 +821,6 @@ def dynamic(self, NoRoutingExecuted): # First, Kinematic/Split routing is solved on all pixels (including MCT pixels) then results are updated # for the MCT pixels. - # Sideflow contribution to MCT grid cells expressed in [m3/s] - SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql - # Grab outflow at the end of the previous routing step t for all pixels) - current state of the MCT pixel ChanQ_0 = self.var.ChanQ.copy() # Outflow (x+dx) at time t (end of previous routing step) (instant) -> used to calc q00 ChanM3_0 = self.var.ChanM3.copy() # Channel storage at time t (end of previous routing step) (instant) V00 @@ -795,6 +831,22 @@ def dynamic(self, NoRoutingExecuted): self.var.ChanM3 = ChanM3 self.var.ChanQAvgDt = ChanQAvgDt # -> used to calc q0m + # # MCT HEADWATER + # self.mctheadwater_module.dynamic_inloop(NoRoutingExecuted) + # # calculate sideflow from MCT headwater pixels + # SideflowChanM3 += self.var.QHeadM3Dt + # # MCT headwater pixels outflow volume per routing sub step [m3] + + # MCT CONFLUENCE + self.mctconfluence_module.dynamic_inloop(NoRoutingExecuted) + # calculate sideflow from MCT confluence pixels + SideflowChanM3 += self.var.QConfM3Dt + # MCT confluence pixels outflow volume per routing sub step [m3] + + # Sideflow contribution to MCT grid cells expressed in [m3/s] + SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql + # SideflowChanMCTM3 = np.where(self.var.IsChannelMCT, SideflowChanM3, 0) + # Solve MCT routing and update current state at MCT pixels self.mct_river_router.routing( ChanQ_0, # -> used to calc q00 From adf248f6736f3ce76897b2adeaab5dbf11305d2e Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 12:49:32 +0200 Subject: [PATCH 03/12] q11 clamp based on physical channel water storage --- src/lisflood/hydrological_modules/mct.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index 7490564c..29e4d424 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/src/lisflood/hydrological_modules/mct.py @@ -336,9 +336,9 @@ def MCTRouting_single( # ql_abs = ql # else: # ql_abs = -ql - # max_q11 = q01 + q00 + ql_abs + V00 / dt - # if q11 > max_q11: - # q11 = max_q11 + max_q11 = q01 + q00 + abs(ql) + V00 / dt + if q11 > max_q11: + q11 = max_q11 # After q11 is computed # max_physical = q01 + q00 + abs(ql) + V00 / dt From 8fd963abb1321f51e1ce57ea7a116746d7cfb676 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 12:51:34 +0200 Subject: [PATCH 04/12] second q11 clamp if deviation from q1mm is too high -> physical q1mm correction --- src/lisflood/hydrological_modules/mct.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index 29e4d424..b48a3a37 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/src/lisflood/hydrological_modules/mct.py @@ -379,12 +379,12 @@ def MCTRouting_single( 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 = 10 * q11 - # elif ratio < 0.1: - # q1mm = 0.1 * 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 From 1c92c960c66cd29fe699891bee84047b03b3846b Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 14:36:06 +0200 Subject: [PATCH 05/12] erging my stability tweaks with cinzias calib points approach to recover the calibration results --- src/lisflood/hydrological_modules/mct.py | 22 ++- .../hydrological_modules/mctconfluence.py | 171 ------------------ src/lisflood/hydrological_modules/routing.py | 40 ++-- 3 files changed, 40 insertions(+), 193 deletions(-) delete mode 100644 src/lisflood/hydrological_modules/mctconfluence.py diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index b48a3a37..8c26185c 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/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. @@ -169,9 +173,19 @@ def mct_routing( q01 = 0.0 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] + # 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) diff --git a/src/lisflood/hydrological_modules/mctconfluence.py b/src/lisflood/hydrological_modules/mctconfluence.py deleted file mode 100644 index fa670825..00000000 --- a/src/lisflood/hydrological_modules/mctconfluence.py +++ /dev/null @@ -1,171 +0,0 @@ -""" - -Copyright 2019 European Union - -Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission -subsequent versions of the EUPL (the "Licence"); - -You may not use this work except in compliance with the Licence. -You may obtain a copy of the Licence at: - -https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt - -Unless required by applicable law or agreed to in writing, -software distributed under the Licence is distributed on an "AS IS" basis, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the Licence for the specific language governing permissions and limitations under the Licence. - -""" - -from __future__ import print_function, absolute_import -from pcraster import Scalar, numpy2pcr, pcr2numpy,downstream, boolean,lddrepair,ifthenelse,lddmask -from nine import range -import numpy as np -from ..global_modules.settings import LisSettings, MaskInfo -from ..global_modules.add1 import loadmap, compressArray, decompress, makenumpy -from ..global_modules.errors import LisfloodWarning -from . import HydroModule - - -class mctconfluence(HydroModule): - """ - Adds contribution from Kinematic cells to MCT cells as a lateral flow when MCT routing is enabled. - This is for MCT cells that have upstream contributions from both kinematic and MCT cells. - - This module handles the initialization and dynamic simulation of Kinematic to MCT cells confluence. - It injects side discharge to the downstream grid cell as lateral flow. - - Attributes: - ----------- - var (object): An object containing all the variables used within the mctconfluence module. - - Methods: - -------- - initial(): Sets up the initial conditions and parameters for the simulation, - including confluence locations. - dynamic_inloop(NoRoutingExecuted: int): Performs dynamic calculations within the routing - loop to simulate the kinematic to MCT confluence. - """ - - module_name = 'MCTConfluence' - - def __init__(self, mctconfluence_variable): - self.var = mctconfluence_variable - - def __init__(self, mctconfluence_variable): - """ - Initializes the MCT confluence module with a given variable object. - - Parameters: - ----------- - mctconfluence_variable: object - An object containing the variables needed for the MCT confluence simulation. - """ - - self.var = mctconfluence_variable - - def initial(self): - """ - Initiates the MCT confluence module by loading the necessary data and maps. - """ - - settings = LisSettings.instance() - option = settings.options - binding = settings.binding - maskinfo = MaskInfo.instance() - if option['MCTRouting']: - - inArPcr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32")) # pcr - # Assign a number to each non-missing pixel as cell id, starting from 0 - inAr = compressArray(inArPcr) #np - - down = (compressArray(downstream(self.var.LddStructuresChan, inArPcr))).astype("int32") # np - # assign to each pixel the cell id of the pixel it is contributing to - # LddStructuresChan do not contain structures, LddStructuresKinematic is same as LddStructuresChan - - LddKinematic = lddmask(self.var.LddStructuresChan, self.var.IsChannelKinematicPcr) - # Mask of LddKinematic with only kinematic cells and no structures - # Sinks are added at the last Kin pixel before confluence with MCT pixels - - maskKinematic = (compressArray(LddKinematic) == 5) #& (self.var.IsUpsOfStructureKinematicC != 1) - - # find location of KIN pixels (only) in LddKin that are upstream of the confluence with an MCT pixel - maskKinematic[compressArray(self.var.AtLastPoint) == 1] = False - # remove sinks that are outlets - # this does NOT include sinks upstream of structures (lakes, reservoirs) and outlets - maskKinematicPcr = boolean(decompress(maskKinematic)) - - self.var.LddChan = ifthenelse(maskKinematicPcr, 5, self.var.LddChan) - self.var.LddKinematic = ifthenelse(maskKinematicPcr, 5, self.var.LddKinematic) - # Adding sinks to Ldd at the last KIN pixels upstream of the confluence with an MCT pixel to LddChan and LddKinematic - # This similar to what is done in structures - # At this point LddChan and LddKinematic have sinks upstteam of structures and of a Kin-MCT confluence and outlets - - self.var.KinematicUpsOfMCTConfluence = np.where(maskKinematic, down, 0) - # find last KIN pixels upstream of the confluence with an MCT pixel and assign it the id of the downstream MCT pixel - - mctconfluence = maskinfo.in_zero() - mctconfluence[np.isin(inAr, self.var.KinematicUpsOfMCTConfluence[self.var.KinematicUpsOfMCTConfluence != 0])] = 1 - # for each element in KinematicUpsOfMCTConfluence (they are Kinematic cells), get the value of the downstream cell, - # then find the position ix of that same value in inAr (vector with numbering of all cells) this is the position of the MCT confluence cell, - # read the inAr value and put 1 in the corrisponding position in mctconfluence - # identify location of MCT pixels that receive a contribution from an upstream KIN pixel - - mctconfluence = np.where(self.var.IsStructureChan, 0, mctconfluence) - # remove the MCT cells with a structure (reservoir/lake) on them - - self.var.MCTConfluenceSitesC = mctconfluence - self.var.MCTConfluenceSitesCC = np.compress(mctconfluence > 0, mctconfluence) - self.var.MCTConfluenceIndex = np.nonzero(mctconfluence)[0] - pass - - - def dynamic_init(self): - """ Initialization of the dynamic part of the MCT confluence module - init mct confluence before sub step routing - """ - settings = LisSettings.instance() - option = settings.options - if option['MCTRouting']: - self.var.QInConfM3Old = np.where(self.var.MCTConfluenceSitesC > 0, self.var.ChanQAvgDt * self.var.DtSec, 0) - - - def dynamic_inloop(self, NoRoutingExecuted: int): - """ - Performs the dynamic simulation of MCT confluence within the routing loop. This method - injects upstream discharge to the downstream grid cell as lateral inflow. - - Parameters: - ----------- - NoRoutingExecuted: integer - The number of routing sub-steps that have been executed. This parameter is used to manage - the accumulation of inflow and outflow over the routing steps. - """ - - settings = LisSettings.instance() - option = settings.options - maskinfo = MaskInfo.instance() - - # self.var.QConfADDEDM3 = maskinfo.in_zero() - - if option['MCTRouting'] and not option['InitLisflood']: - - InvDtSecDay = 1 / float(86400) - # InvDtSecDay=self.var.InvDtSec - - lateralflow = np.bincount(self.var.KinematicUpsOfMCTConfluence, weights=self.var.ChanQAvgDt)[self.var.MCTConfluenceIndex] #same as Qin - # contribution to the MCT pixel from upstream Kinematic pixels - - self.var.QInConfM3 = maskinfo.in_zero() - np.put(self.var.QInConfM3, self.var.MCTConfluenceIndex, lateralflow * self.var.DtSec) - self.var.QDeltaConfM3 = (self.var.QInConfM3 - self.var.QInConfM3Old) * self.var.InvNoRoutSteps - # difference between old and new lateral flow per sub step - # in order to calculate the amount of lateral flow in the routing loop - - self.var.QConfM3Dt = (self.var.QInConfM3Old + (NoRoutingExecuted + 1) * self.var.QDeltaConfM3) * self.var.InvNoRoutSteps - # output to the MCT confluence cell - - self.var.QInConfM3Old = self.var.QInConfM3.copy() - # save the lateral flow for next step - - # self.var.QConfADDEDM3 += self.var.QConfM3Dt \ No newline at end of file diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index 2adc1074..14012a63 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -16,7 +16,7 @@ """ from __future__ import print_function, absolute_import -from pcraster import lddmask, accuflux, boolean, scalar, downstream, pit, path, lddrepair, ifthenelse, cover, nominal, uniqueid, \ +from pcraster import lddmask, accuflux, boolean, downstream, pit, path, lddrepair, ifthenelse, cover, nominal, uniqueid, \ catchment, upstream, pcr2numpy import warnings @@ -30,7 +30,6 @@ from .transmission import transmission from .kinematic_wave_parallel import kinematicWave, kwpt from .mct import MCTWave -from .mctconfluence import mctconfluence from ..global_modules.settings import LisSettings, MaskInfo from ..global_modules.errors import LisfloodWarning @@ -51,7 +50,8 @@ class routing(HydroModule): 'ChanBottomWMult', 'ChanDepthTMult', 'ChanSMult'], 'SplitRouting': ['CrossSection2AreaInitValue', 'PrevSideflowInitValue', 'CalChanMan2'], 'dynamicWave': ['ChannelsDynamic'], - 'MCTRouting': ['ChannelsMCT', 'ChanGradMaxMCT', 'PrevCmMCTInitValue', 'PrevDmMCTInitValue', 'CalChanMan3']} + 'MCTRouting': ['ChannelsMCT', 'ChanGradMaxMCT', 'PrevCmMCTInitValue', 'PrevDmMCTInitValue', 'CalChanMan3'], + 'simulateCalibrationPoints': ['CalibrationPoints']} module_name = 'Routing' def __init__(self, routing_variable): @@ -62,7 +62,6 @@ def __init__(self, routing_variable): self.polder_module = polder(self.var) self.inflow_module = inflow(self.var) self.transmission_module = transmission(self.var) - self.mctconfluence_module = mctconfluence(self.var) # -------------------------------------------------------------------------- # -------------------------------------------------------------------------- @@ -291,7 +290,7 @@ def initial(self): # set BankFullPerc to 0.5 for half bankfull # Channel volume initialization for MCT cells - TotalCrossSectionAreaHalfBankFull = np.where(self.var.IsChannelKinematic, TotalCrossSectionAreaHalfBankFull, 0.01 * self.var.TotalCrossSectionAreaBankFull) + # TotalCrossSectionAreaHalfBankFull = np.where(self.var.IsChannelKinematic, TotalCrossSectionAreaHalfBankFull, 0.01 * self.var.TotalCrossSectionAreaBankFull) # set initial volume in MCT cells to 1% of bankfull TotalCrossSectionAreaInitValue = loadmap('TotalCrossSectionAreaInitValue') @@ -641,6 +640,19 @@ def initialMCT(self): self.var.PrevDm0 = np.where(PrevDmMCT == -9999, maskinfo.in_zero(), PrevDmMCT) #np # Reynolds number (Dm) for MCT at previous time step t0 + # ************************************************************ + # ***** CALIBRATION POINTS ******** + # ************************************************************ + CalibPoints = maskinfo.in_zero() + if option['simulateCalibrationPoints']: + CalibPoints = loadmap('CalibrationPoints') # 1D array size all catchment pixels + # read location of calibration points + + inAr = np.arange(maskinfo.info.mapC[0], dtype="int32") # np + # Assign a number to each non-missing pixel as cell id, by row starting from 0 + CalibPointsIds = inAr[CalibPoints > 0] + # pixel id of calibration points + # ************************************************************ # ***** INITIALISE MUSKINGUM-CUNGE-TODINI WAVE ROUTER ******** @@ -648,6 +660,9 @@ def initialMCT(self): mct_ldd = self.compress_mct(compressArray(self.var.LddMCT)) # Compress LddMCT to array with MCT pixels only + # mct_CalInflowPoints = self.compress_mct(self.var.CalInflowPoints) + # # Compress CalInflowPoints to array with MCT pixels only + mapping_mct = self.compress_mct(range(len(self.var.ChanLength))) # create mapping from global domain pixels index to MCT pixels index @@ -661,7 +676,8 @@ def initialMCT(self): self.var.ChanSdXdY, # Riverbed side slope self.var.DtRouting, # computation time step for routing [s] self.river_router, # class - mapping_mct # MCT pixels mapping + mapping_mct, # MCT pixels mapping + CalibPointsIds, # id of calibrationn points in full LDD ) @@ -831,18 +847,6 @@ def dynamic(self, NoRoutingExecuted): self.var.ChanM3 = ChanM3 self.var.ChanQAvgDt = ChanQAvgDt # -> used to calc q0m - # # MCT HEADWATER - # self.mctheadwater_module.dynamic_inloop(NoRoutingExecuted) - # # calculate sideflow from MCT headwater pixels - # SideflowChanM3 += self.var.QHeadM3Dt - # # MCT headwater pixels outflow volume per routing sub step [m3] - - # MCT CONFLUENCE - self.mctconfluence_module.dynamic_inloop(NoRoutingExecuted) - # calculate sideflow from MCT confluence pixels - SideflowChanM3 += self.var.QConfM3Dt - # MCT confluence pixels outflow volume per routing sub step [m3] - # Sideflow contribution to MCT grid cells expressed in [m3/s] SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql # SideflowChanMCTM3 = np.where(self.var.IsChannelMCT, SideflowChanM3, 0) From 3e658fb9ce4700c37567dd0d7167691f3b412574 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 14:41:03 +0200 Subject: [PATCH 06/12] erging my stability tweaks with cinzias calib points approach to recover the calibration results II --- src/lisflood/Lisflood_initial.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index 084377dd..656e2e6b 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/src/lisflood/Lisflood_initial.py @@ -218,11 +218,6 @@ def __init__(self): # 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.mctheadwater_module.initial() - # # initialising MCT headwater pixels and adding pixels upstream of MCT headwater to the LDD - self.mctconfluence_module.initial() - # initialising MCT confluence points and adding MCT confluence sinks to the LDD - # ---------------------------------------------------------------------- # ---------------------------------------------------------------------- @@ -238,9 +233,6 @@ def __init__(self): self.routing_module.initialMCT() # initialising Muskingum-Cunge-Todini routing for channel - # self.mctheadwater_module.dynamic_init() - self.mctconfluence_module.dynamic_init() - # self.routing_module.initialKinematicWave() # this cannot be here because I need river_router in the MCT initialization From f90a3ad6f77ff378bae147c2d89549c51d2bdfce Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 14:42:11 +0200 Subject: [PATCH 07/12] erging my stability tweaks with cinzias calib points approach to recover the calibration results III --- src/lisflood/Lisflood_initial.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index 656e2e6b..a1f29eb4 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/src/lisflood/Lisflood_initial.py @@ -139,8 +139,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.mctheadwater_module = mctheadwater(self) - self.mctconfluence_module = mctconfluence(self) + self.polder_module = polder(self) self.waterabstraction_module = waterabstraction(self) self.indicatorcalc_module = indicatorcalc(self) From 2a892c863b52f4b6698e030704137c7ebf82e912 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 14:43:17 +0200 Subject: [PATCH 08/12] erging my stability tweaks with cinzias calib points approach to recover the calibration results IV --- src/lisflood/Lisflood_initial.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index a1f29eb4..caa48104 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/src/lisflood/Lisflood_initial.py @@ -49,7 +49,6 @@ from .hydrological_modules.surface_routing import surface_routing from .hydrological_modules.reservoir import Reservoir from .hydrological_modules.lakes import lakes -from .hydrological_modules.mctconfluence import mctconfluence from .hydrological_modules.polder import polder from .hydrological_modules.waterabstraction import waterabstraction from .hydrological_modules.indicatorcalc import indicatorcalc From da7bf793fe9938252090d3d9634253a6a0fddb36 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 15:00:58 +0200 Subject: [PATCH 09/12] erging also the simulateCalibrationPoints in default_options.py --- src/lisflood/global_modules/default_options.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py index 94181f80..aa2a9968 100644 --- a/src/lisflood/global_modules/default_options.py +++ b/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, From 88885930f8b388f862c7ec165d9e15d3f062abcb Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Wed, 20 May 2026 18:49:01 +0200 Subject: [PATCH 10/12] adding one more change from cinzia calib points --- src/lisflood/hydrological_modules/mct.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index 8c26185c..15c14fed 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/src/lisflood/hydrological_modules/mct.py @@ -171,13 +171,14 @@ 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 ##################################################################################################### # 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] + ql += ChanQAvgDt[ups_pix] # avoid += # Sideflow during step dt including contribution from calibration pixel else: # not a calibration point - go as usual @@ -195,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 From 3dbd73503461ab3ec4876b99313ae75e2544fc77 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Fri, 22 May 2026 11:45:29 +0200 Subject: [PATCH 11/12] remove q1mm aggressiv clamp --- src/lisflood/hydrological_modules/mct.py | 18 ++++++------------ src/lisflood/hydrological_modules/routing.py | 4 ++-- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py index 15c14fed..2a39f339 100644 --- a/src/lisflood/hydrological_modules/mct.py +++ b/src/lisflood/hydrological_modules/mct.py @@ -355,12 +355,6 @@ def MCTRouting_single( if q11 > max_q11: q11 = max_q11 - # After q11 is computed - # max_physical = q01 + q00 + abs(ql) + V00 / dt - # if q11 > max_physical * 1.1: # exceeds physical bound by >10% - # print("MASS VIOLATION q11=", q11, "max_phys=", max_physical, - # "q01=", q01, "q00=", q00, "ql=", ql, "V00=", V00) - #### end of for loop # # cmcheck @@ -394,12 +388,12 @@ def MCTRouting_single( 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 + # 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 diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index 14012a63..886ef50c 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -868,7 +868,7 @@ def dynamic(self, NoRoutingExecuted): # Case 1: ChanQAvgDt too small relative to ChanQ (flood arrival) mct_avg_too_small = ( self.var.IsChannelMCT & - (self.var.ChanQ > 100.) & + (self.var.ChanQ > 0.1) & (self.var.ChanQAvgDt < 0.1 * self.var.ChanQ) ) self.var.ChanQAvgDt = np.where(mct_avg_too_small, self.var.ChanQ, self.var.ChanQAvgDt) @@ -876,7 +876,7 @@ def dynamic(self, NoRoutingExecuted): # Case 2: ChanQAvgDt too large relative to ChanQ (flood recession / cold start) mct_avg_too_large = ( self.var.IsChannelMCT & - (self.var.ChanQAvgDt > 100.) & + (self.var.ChanQAvgDt > 0.1) & (self.var.ChanQ < 0.1 * self.var.ChanQAvgDt) ) self.var.ChanQAvgDt = np.where(mct_avg_too_large, self.var.ChanQ, self.var.ChanQAvgDt) From b155c8c65f229fc560c1279d334418bd5308b5ad Mon Sep 17 00:00:00 2001 From: doc78 Date: Fri, 29 May 2026 08:04:53 +0000 Subject: [PATCH 12/12] Added warmstart UTs with calibpoints for ETRS89. Updated temp version number --- VERSION | 2 +- .../LF_ETRS89_UseCase/settings/mct_cold.xml | 15 +++++ .../LF_ETRS89_UseCase/settings/mct_warm.xml | 15 +++++ tests/test_mct_warmstart_slow.py | 55 +++++++++++++++++++ 4 files changed, 86 insertions(+), 1 deletion(-) diff --git a/VERSION b/VERSION index c3ef9b53..9286fc2e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -4.3.1.505 \ No newline at end of file +4.3.1.506 \ No newline at end of file diff --git a/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml b/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml index e4b4879f..91caf0ce 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/mct_cold.xml @@ -47,6 +47,7 @@ You can use builtin path variables in this template and reference to other paths + @@ -1681,6 +1682,13 @@ You can use builtin path variables in this template and reference to other paths + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + + OPTIONAL: observed or simulated input hydrographs as @@ -5476,6 +5484,13 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + + Observed or simulated input hydrographs as diff --git a/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml b/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml index 639730fc..90d26839 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/mct_warm.xml @@ -47,6 +47,7 @@ You can use builtin path variables in this template and reference to other paths + @@ -1650,6 +1651,13 @@ You can use builtin path variables in this template and reference to other paths + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + + OPTIONAL: observed or simulated input hydrographs as @@ -5426,6 +5434,13 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT inflow hydrographs + + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + diff --git a/tests/test_mct_warmstart_slow.py b/tests/test_mct_warmstart_slow.py index 88d26fd5..baa6b3c8 100644 --- a/tests/test_mct_warmstart_slow.py +++ b/tests/test_mct_warmstart_slow.py @@ -77,6 +77,42 @@ def test_mct_warmstart_6h(self): report_steps = '38220..38830' self.run_warmstart_by_dtsec('mct_all', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps) + def test_mct_only_warmstart_daily_withcalibpoints(self): + calendar_day_start = '02/01/1990 06:00' + step_start = '02/01/2016 06:00' + step_end = '31/12/2016 06:00' + dt_sec = 86400 + dt_sec_channel = 3600 + report_steps = '9496..9861' + self.run_warmstart_by_dtsec('mct_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps) + + def test_mct_warmstart_daily_withcalibpoints(self): + calendar_day_start = '02/01/1990 06:00' + step_start = '02/01/2016 06:00' + step_end = '31/12/2016 06:00' + dt_sec = 86400 + dt_sec_channel = 3600 + report_steps = '9496..9861' + self.run_warmstart_by_dtsec('mct_all_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps) + + def test_mct_only_warmstart_6h_withcalibpoints(self): + calendar_day_start = '02/01/1990 06:00' + step_start = '01/03/2016 06:00' + step_end = '31/07/2016 06:00' + dt_sec = 21600 + dt_sec_channel = 3600 + report_steps = '38220..38830' + self.run_warmstart_by_dtsec('mct_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps) + + def test_mct_warmstart_6h_withcalibpoints(self): + calendar_day_start = '02/01/1990 06:00' + step_start = '01/03/2016 06:00' + step_end = '31/07/2016 06:00' + dt_sec = 21600 + dt_sec_channel = 3600 + report_steps = '38220..38830' + self.run_warmstart_by_dtsec('mct_all_calibpoints', dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps=report_steps) + def run_warmstart_by_dtsec(self, mct_case, dt_sec, dt_sec_channel, step_end, step_start, calendar_day_start,report_steps='1..9999'): mk_path_out(os.path.join(self.case_dir, 'out')) @@ -101,6 +137,25 @@ def run_warmstart_by_dtsec(self, mct_case, dt_sec, dt_sec_channel, step_end, ste ] opts_to_unset = ['repMBTs', 'simulateReservoirs'] + elif mct_case == 'mct_calibpoints': + opts_to_set = ['repStateMaps','repDischargeMaps', + 'TransLoss', + 'simulateCalibrationPoints'] + opts_to_unset = ['repMBTs', 'simulateReservoirs', 'simulateLakes'] + elif mct_case == 'mct_all_calibpoints': + opts_to_set = ['repStateMaps', + 'repDischargeMaps', + 'wateruse', + 'drainedIrrigation', + 'riceIrrigation', + 'openwaterevapo', + 'simulateLakes', + 'TransLoss', + 'simulateCalibrationPoints' + ] + opts_to_unset = ['repMBTs', + 'simulateReservoirs'] + settings_longrun = setoptions(self.settings_files['cold'], opts_to_set=opts_to_set,