Refactor routing

docs/changelog.qmd

@@ -16,6 +16,11 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 - Mutating BMI model control functions (`update`, `update_until` and `finalize`) and
   extended mutating BMI functions (`load_state` and `save_state`) should return `nothing`.
 - Added downloading of testdata to Dockerfile, to ensure an image was able to build.
+- The reservoir (`reservoir_index_f`) and lake (`lake_index_f`) indices as part of
+  `network.river` were not correct. These were mapped to their own index in the
+  `SimpleReservoir` and `Lake` struct, and not to the corresponding river index. This
+  resulted in incorrect surface water abstractions from reservoir and lake volumes, and
+  surface water abstractions were set at zero at the wrong river locations.
 
 ### Changed
 - Removed vertical concepts `HBV` and `FLEXTopo`.
@@ -31,6 +36,13 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   model `variables`, `parameters` and `boundary_conditions` (if applicable), including
   associated functions for initializing and updating these model components. The original
   long update function of the `SBM` soil part has been split into separate functions.
+- Refactor the lateral (routing) components: as for the vertical `SBM` concept split the
+  structs into `variables`, `parameters` and `boundary_conditions` (if applicable).
+- Timestepping method parameters for solving the kinematic wave and local inertial
+  approaches for river and overland flow are moved to a `TimeStepping` struct. The
+  timestepping implementation for the kinematic wave is now similar to the local inertial
+  method: a stable timestep is computed for each sub timestep (or a fixed sub timestep is
+  used) as part of a while loop (for each model timestep).
 
 ### Added
 - Support direct output of snow and glacier melt, and add computation of snow water

server/test/client.jl

@@ -32,13 +32,13 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 203)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 203)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 202)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 202)
     to_check = [
         "vertical.soil.parameters.nlayers",
         "vertical.soil.parameters.theta_r",
-        "lateral.river.q",
-        "lateral.river.reservoir.outflow",
+        "lateral.river.variables.q",
+        "lateral.river.boundary_conditions.reservoir.variables.outflow",
     ]
     retrieved_vars = request((fn = "get_input_var_names",))["input_var_names"]
     @test all(x -> x in retrieved_vars, to_check)
@@ -49,12 +49,12 @@ end
 zi_size = 0
 vwc_1_size = 0
 @testset "variable information and get and set functions" begin
-    @test request((fn = "get_var_itemsize", name = "lateral.subsurface.ssf")) ==
+    @test request((fn = "get_var_itemsize", name = "lateral.subsurface.variables.ssf")) ==
           Dict("var_itemsize" => sizeof(Wflow.Float))
     @test request((fn = "get_var_type", name = "vertical.n"))["status"] == "ERROR"
     @test request((fn = "get_var_units", name = "vertical.soil.parameters.theta_s")) ==
           Dict("var_units" => "-")
-    @test request((fn = "get_var_location", name = "lateral.river.q")) ==
+    @test request((fn = "get_var_location", name = "lateral.river.variables.q")) ==
           Dict("var_location" => "node")
     zi_nbytes =
         request((fn = "get_var_nbytes", name = "vertical.soil.variables.zi"))["var_nbytes"]
@@ -68,19 +68,20 @@ vwc_1_size = 0
     vwc_1_itemsize =
         request((fn = "get_var_itemsize", name = "vertical.soil.variables.vwc[1]"))["var_itemsize"]
     vwc_1_size = Int(vwc_1_nbytes / vwc_1_itemsize)
-    @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
+    @test request((fn = "get_var_grid", name = "lateral.river.variables.h")) ==
+          Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.soil.variables.zi", dest = fill(0.0, zi_size))
     @test mean(request(msg)["value"]) ≈ 277.3620724821974
     msg = (fn = "get_value_ptr", name = "vertical.soil.parameters.theta_s")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
         fn = "get_value_at_indices",
-        name = "lateral.river.q",
+        name = "lateral.river.variables.q",
         dest = [0.0, 0.0, 0.0],
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.198747900215207f0, 2.6880427720508515f0, 3.4848783702629564f0]
+          [2.1007361866518766, 2.5702292750107687, 3.2904803551115727]
     msg =
         (fn = "set_value", name = "vertical.soil.variables.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")

server/test/sbm_config.toml

@@ -31,18 +31,18 @@ ustorelayerdepth = "ustorelayerdepth"
 snow_storage = "snow"
 snow_water = "snowwater"
 
-[state.lateral.river]
+[state.lateral.river.variables]
 h = "h_river"
 h_av = "h_av_river"
 q = "q_river"
 
-[state.lateral.river.reservoir]
+[state.lateral.river.boundary_conditions.reservoir.variables]
 volume = "volume_reservoir"
 
-[state.lateral.subsurface]
+[state.lateral.subsurface.variables]
 ssf = "ssf"
 
-[state.lateral.land]
+[state.lateral.land.variables]
 h = "h_land"
 h_av = "h_av_land"
 q = "q_land"
@@ -112,7 +112,7 @@ offset = 0.0
 
 [input.lateral.river]
 length = "wflow_riverlength"
-n = "N_River"
+mannings_n = "N_River"
 slope = "RiverSlope"
 width = "wflow_riverwidth"
 bankfull_elevation = "RiverZ"
@@ -132,7 +132,7 @@ targetminfrac = "ResTargetMinFrac"
 ksathorfrac = "KsatHorFrac"
 
 [input.lateral.land]
-n = "N"
+mannings_n = "N"
 slope = "Slope"
 
 [model]
@@ -161,17 +161,17 @@ ustorelayerdepth = "ustorelayerdepth"
 snow_storage = "snow"
 snow_water = "snowwater"
 
-[output.lateral.river]
+[output.lateral.river.variables]
 h = "h_river"
 q = "q_river"
 
-[output.lateral.river.reservoir]
+[output.lateral.river.boundary_conditions.reservoir.variables]
 volume = "volume_reservoir"
 
-[output.lateral.subsurface]
+[output.lateral.subsurface.variables]
 ssf = "ssf"
 
-[output.lateral.land]
+[output.lateral.land.variables]
 h = "h_land"
 q = "q_land"
 
@@ -181,7 +181,7 @@ path = "output_scalar_moselle.nc"
 [[netcdf.variable]]
 name = "Q"
 map = "gauges"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 
 [[netcdf.variable]]
 coordinate.x = 6.255
@@ -202,13 +202,13 @@ path = "output_moselle.csv"
 
 [[csv.column]]
 header = "Q"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 reducer = "maximum"
 
 [[csv.column]]
 header = "volume"
 index = 1
-parameter = "lateral.river.reservoir.volume"
+parameter = "lateral.river.boundary_conditions.reservoir.variables.volume"
 
 [[csv.column]]
 coordinate.x = 6.255
@@ -232,7 +232,7 @@ parameter = "vertical.atmospheric_forcing.temperature"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.q"
+parameter = "lateral.river.variables.q"
 
 [[csv.column]]
 header = "recharge"
@@ -255,8 +255,19 @@ components = [
   "vertical.snow.boundary_conditions",
   "vertical.snow.variables",
   "vertical.snow.parameters",
-  "lateral.subsurface",
-  "lateral.land",
-  "lateral.river",
-  "lateral.river.reservoir",
+  "lateral.subsurface.boundary_conditions",
+  "lateral.subsurface.variables",
+  "lateral.subsurface.parameters",
+  "lateral.subsurface.parameters.kh_profile",
+  "lateral.land.boundary_conditions",
+  "lateral.land.variables",
+  "lateral.land.variables.flow",
+  "lateral.land.parameters",
+  "lateral.river.variables",
+  "lateral.river.parameters",
+  "lateral.river.parameters.flow",
+  "lateral.river.boundary_conditions",
+  "lateral.river.boundary_conditions.reservoir.boundary_conditions",
+  "lateral.river.boundary_conditions.reservoir.parameters",
+  "lateral.river.boundary_conditions.reservoir.variables",
 ]

src/Wflow.jl

@@ -49,7 +49,7 @@ using Parameters: @with_kw
 using Polyester: @batch
 using ProgressLogging: @progress
 using StaticArrays: SVector, pushfirst, setindex
-using Statistics: mean, median, quantile!
+using Statistics: mean, median, quantile!, quantile
 using TerminalLoggers
 using TOML: TOML
 
@@ -141,8 +141,17 @@ Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))
 
 include("forcing.jl")
 include("parameters.jl")
-include("flow.jl")
-include("horizontal_process.jl")
+include("groundwater/connectivity.jl")
+include("groundwater/aquifer.jl")
+include("groundwater/boundary_conditions.jl")
+include("routing/timestepping.jl")
+include("routing/subsurface.jl")
+include("routing/reservoir.jl")
+include("routing/lake.jl")
+include("routing/surface_kinwave.jl")
+include("routing/surface_local_inertial.jl")
+include("routing/surface_routing.jl")
+include("routing/routing_process.jl")
 include("vegetation/rainfall_interception.jl")
 include("vegetation/canopy.jl")
 include("snow/snow_process.jl")
@@ -155,12 +164,8 @@ include("soil/soil_process.jl")
 include("sbm.jl")
 include("demand/water_demand.jl")
 include("sediment.jl")
-include("reservoir_lake.jl")
 include("sbm_model.jl")
 include("sediment_model.jl")
-include("groundwater/connectivity.jl")
-include("groundwater/aquifer.jl")
-include("groundwater/boundary_conditions.jl")
 include("sbm_gwf_model.jl")
 include("utils.jl")
 include("bmi.jl")

src/bmi.jl

@@ -153,7 +153,7 @@ function BMI.get_var_grid(model::Model, name::String)
     s = split(name, "[")
     key = symbols(first(s))
     if exchange(param(model, key))
-        type = typeof(param(model, key[1:(end - 1)]))
+        type = typeof(param(model, key[1:2]))
         return if :reservoir in key
             0
         elseif :lake in key
@@ -162,9 +162,9 @@ function BMI.get_var_grid(model::Model, name::String)
             2
         elseif :river in key
             3
-        elseif type <: ShallowWaterLand && occursin("x", s[end])
+        elseif type <: LocalInertialOverlandFlow && occursin("x", s[end])
             4
-        elseif type <: ShallowWaterLand && occursin("y", s[end])
+        elseif type <: LocalInertialOverlandFlow && occursin("y", s[end])
             5
         else
             6
@@ -374,7 +374,7 @@ function BMI.get_grid_edge_nodes(model::Model, grid::Int, edge_nodes::Vector{Int
     m = div(n, 2)
     # inactive nodes (boundary/ghost points) are set at -999
     if grid == 3
-        nodes_at_edge = adjacent_nodes_at_link(network.river.graph)
+        nodes_at_edge = adjacent_nodes_at_edge(network.river.graph)
         nodes_at_edge.dst[nodes_at_edge.dst .== m + 1] .= -999
         edge_nodes[range(1, n; step = 2)] = nodes_at_edge.src
         edge_nodes[range(2, n; step = 2)] = nodes_at_edge.dst