From 82b72797e1f4dc4960e6a072e1253b3c2b8d1b51 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Wed, 13 Sep 2023 16:37:11 +0200 Subject: [PATCH] Only gravity driven `Outlet` flow (#549) Fixes https://github.com/Deltares/Ribasim/issues/548. Currently the outlet switches on and off in a discrete way due to a switch in sign of the level difference. This can result in the outlet rapidly turning on and off: ![flipperende_outlet](https://github.com/Deltares/Ribasim/assets/74617371/d834ffa9-ac48-4fb1-99c5-732369f43c35) To solve this problem, a reduction factor should be applied to the outlet flow, based on the level difference. This reduction factor will have to be taken into account in the implementation of the analytical Jacobian (with an unsatisfying effort to improvement ratio). This I think will be a recurring problem until we are able to switch to an AD Jacobian. --------- Co-authored-by: Martijn Visser --- build/libribasim/src/libribasim.jl | 6 +- core/src/bmi.jl | 6 +- core/src/config.jl | 8 +- core/src/create.jl | 8 +- core/src/io.jl | 6 +- core/src/solve.jl | 121 +++++++++--- core/src/utils.jl | 35 +++- core/src/validation.jl | 1 + core/test/control.jl | 4 +- core/test/libribasim.jl | 6 +- core/test/run_models.jl | 29 ++- core/test/utils.jl | 10 + core/test/validation.jl | 1 + docs/core/equations.qmd | 175 ++++++++---------- docs/core/usage.qmd | 22 ++- docs/python/examples.ipynb | 6 +- docs/schema/OutletStatic.schema.json | 7 + python/ribasim/ribasim/geometry/node.py | 2 +- python/ribasim/ribasim/model.py | 1 + python/ribasim/ribasim/models.py | 3 +- .../ribasim_testmodels/__init__.py | 2 + .../ribasim_testmodels/basic.py | 107 ++++++++++- .../ribasim_testmodels/pid_control.py | 4 +- qgis/core/nodes.py | 6 +- 24 files changed, 405 insertions(+), 171 deletions(-) diff --git a/build/libribasim/src/libribasim.jl b/build/libribasim/src/libribasim.jl index d4551e159..ccaf3fb54 100644 --- a/build/libribasim/src/libribasim.jl +++ b/build/libribasim/src/libribasim.jl @@ -27,7 +27,7 @@ This expands to ``` try global model - isnothing(model) && error("Model not initialized") + model === nothing && error("Model not initialized") BMI.update(model) catch Base.invokelatest(Base.display_error, current_exceptions()) @@ -40,7 +40,7 @@ macro try_c(ex) return quote try global model - isnothing(model) && error("Model not initialized") + model === nothing && error("Model not initialized") $(esc(ex)) catch e global last_error_message @@ -81,7 +81,7 @@ end Base.@ccallable function finalize()::Cint @try_c_uninitialized begin - if !isnothing(model) + if model !== nothing BMI.finalize(model) end model = nothing diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 2890ead2c..7aa8b5497 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -53,7 +53,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model for id in pid_control.node_id id_controlled = only(outneighbors(connectivity.graph_control, id)) pump_idx = findsorted(pump.node_id, id_controlled) - if isnothing(pump_idx) + if pump_idx === nothing outlet_idx = findsorted(outlet.node_id, id_controlled) outlet.is_pid_controlled[outlet_idx] = true else @@ -247,7 +247,7 @@ function get_value( if hasindex_basin _, level = get_area_and_level(basin, basin_idx, u[basin_idx]) - elseif !isnothing(level_boundary_idx) + elseif level_boundary_idx !== nothing level = level_boundary.level[level_boundary_idx](t + Δt) else error( @@ -260,7 +260,7 @@ function get_value( elseif variable == "flow_rate" flow_boundary_idx = findsorted(flow_boundary.node_id, feature_id) - if isnothing(flow_boundary_idx) + if flow_boundary_idx === nothing error("Flow condition node #$feature_id is not a flow boundary.") end diff --git a/core/src/config.jl b/core/src/config.jl index fee738f45..6073a2187 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -84,7 +84,7 @@ end function Base.convert(::Type{Compression}, str::AbstractString) i = findfirst(==(Symbol(str)) ∘ Symbol, instances(Compression)) - if isnothing(i) + if i === nothing throw( ArgumentError( "Compression algorithm $str not supported, choose one of: $(join(instances(Compression), " ")).", @@ -149,7 +149,7 @@ function Base.show(io::IO, c::Config) println(io, "Ribasim Config") for field in fieldnames(typeof(c)) f = getfield(c, field) - isnothing(f) || println(io, "\t$field\t= $f") + f === nothing || println(io, "\t$field\t= $f") end end @@ -157,7 +157,7 @@ function Base.show(io::IO, c::TableOption) first = true for field in fieldnames(typeof(c)) f = getfield(c, field) - if !isnothing(f) + if f !== nothing first && (first = false; println(io)) println(io, "\t\t$field\t= $f") end @@ -180,7 +180,7 @@ const algorithms = Dict{String, Type}( "Create an OrdinaryDiffEqAlgorithm from solver config" function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm algotype = get(algorithms, solver.algorithm, nothing) - if isnothing(algotype) + if algotype === nothing options = join(keys(algorithms), ", ") error("Given solver algorithm $(solver.algorithm) not supported.\n\ Available options are: ($(options)).") diff --git a/core/src/create.jl b/core/src/create.jl index c75a6c0e8..7988fb6f3 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -74,14 +74,14 @@ function parse_static_and_time( end # Get node IDs of static nodes if the static table exists - static_node_ids = if isnothing(static) + static_node_ids = if static === nothing Set{Int}() else Set(static.node_id) end # Get node IDs of transient nodes if the time table exists - time_node_ids = if isnothing(time) + time_node_ids = if time === nothing Set{Int}() else Set(time.node_id) @@ -433,7 +433,8 @@ end function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet static = load_structvector(db, config, OutletStaticV1) - defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true) + defaults = + (; min_flow_rate = 0.0, max_flow_rate = NaN, min_crest_level = NaN, active = true) parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults) is_pid_controlled = falses(length(parsed_parameters.node_id)) @@ -454,6 +455,7 @@ function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet flow_rate, parsed_parameters.min_flow_rate, parsed_parameters.max_flow_rate, + parsed_parameters.min_crest_level, parsed_parameters.control_mapping, is_pid_controlled, ) diff --git a/core/src/io.jl b/core/src/io.jl index 4be67f2db..afe50db34 100644 --- a/core/src/io.jl +++ b/core/src/io.jl @@ -50,7 +50,7 @@ function load_data( path = getfield(getfield(config, snake_case(node)), kind) sqltable = tablename(schema) - table = if !isnothing(path) + table = if path !== nothing table_path = input_path(config, path) Table(read(table_path)) elseif exists(db, sqltable) @@ -76,7 +76,7 @@ function load_structvector( )::StructVector{T} where {T <: AbstractRow} table = load_data(db, config, T) - if isnothing(table) + if table === nothing return StructVector{T}(undef, 0) end @@ -98,7 +98,7 @@ function load_structvector( table = StructVector{T}(nt) sv = Legolas._schema_version_from_record_type(T) tableschema = Tables.schema(table) - if declared(sv) && !isnothing(tableschema) + if declared(sv) && tableschema !== nothing validate(tableschema, sv) # R = Legolas.record_type(sv) # foreach(R, Tables.rows(table)) # construct each row diff --git a/core/src/solve.jl b/core/src/solve.jl index 6e28a7f32..634acdab0 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -318,6 +318,7 @@ struct Outlet{T} <: AbstractParameterNode flow_rate::T min_flow_rate::Vector{Float64} max_flow_rate::Vector{Float64} + min_crest_level::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} is_pid_controlled::BitVector @@ -327,6 +328,7 @@ struct Outlet{T} <: AbstractParameterNode flow_rate::T, min_flow_rate, max_flow_rate, + min_crest_level, control_mapping, is_pid_controlled, ) where {T} @@ -337,6 +339,7 @@ struct Outlet{T} <: AbstractParameterNode flow_rate, min_flow_rate, max_flow_rate, + min_crest_level, control_mapping, is_pid_controlled, ) @@ -501,7 +504,7 @@ function set_current_basin_properties!( end """ -Linearize the evaporation flux when at small water depths +Smoothly let the evaporation flux go to 0 when at small water depths Currently at less than 0.1 m. """ function formulate!( @@ -520,12 +523,12 @@ function formulate!( bottom = basin.level[i][1] fixed_area = basin.area[i][end] depth = max(level - bottom, 0.0) - reduction_factor = clamp(depth, 0.0, 0.1) / 0.1 + factor = reduction_factor(depth, 0.1) precipitation = fixed_area * basin.precipitation[i] - evaporation = area * reduction_factor * basin.potential_evaporation[i] + evaporation = area * factor * basin.potential_evaporation[i] drainage = basin.drainage[i] - infiltration = reduction_factor * basin.infiltration[i] + infiltration = factor * basin.infiltration[i] du.storage[i] += precipitation - evaporation + drainage - infiltration end @@ -579,7 +582,7 @@ function continuous_control!( if !active[i] du.integral[i] = 0.0 u.integral[i] = 0.0 - return + continue end du.integral[i] = pid_error[i] @@ -590,11 +593,29 @@ function continuous_control!( controlled_node_id = only(outneighbors(graph_control, id)) controls_pump = (controlled_node_id in pump.node_id) + # No flow of outlet if source level is lower than target level + if !controls_pump + src_id = only(inneighbors(graph_flow, controlled_node_id)) + dst_id = only(outneighbors(graph_flow, controlled_node_id)) + + src_level = get_level(p, src_id, current_level, t) + dst_level = get_level(p, dst_id, current_level, t) + + if src_level === nothing || dst_level === nothing + factor_outlet = 1.0 + else + Δlevel = src_level - dst_level + factor_outlet = reduction_factor(Δlevel, 0.1) + end + else + factor_outlet = 1.0 + end + if controls_pump controlled_node_idx = findsorted(pump.node_id, controlled_node_id) listened_basin_storage = u.storage[listened_node_idx] - reduction_factor = clamp(listened_basin_storage, 0.0, 10.0) / 10.0 + factor_basin = reduction_factor(listened_basin_storage, 10.0) else controlled_node_idx = findsorted(outlet.node_id, controlled_node_id) @@ -603,12 +624,13 @@ function continuous_control!( has_index, upstream_basin_idx = id_index(basin.node_id, upstream_node_id) if has_index upstream_basin_storage = u.storage[upstream_basin_idx] - reduction_factor = clamp(upstream_basin_storage, 0.0, 10.0) / 10.0 + factor_basin = reduction_factor(upstream_basin_storage, 10.0) else - reduction_factor = 1.0 + factor_basin = 1.0 end end + factor = factor_basin * factor_outlet flow_rate = 0.0 K_p, K_i, K_d = pid_params[i](t) @@ -616,17 +638,17 @@ function continuous_control!( if !iszero(K_d) # dlevel/dstorage = 1/area area = current_area[listened_node_idx] - D = 1.0 - K_d * reduction_factor / area + D = 1.0 - K_d * factor / area else D = 1.0 end if !iszero(K_p) - flow_rate += reduction_factor * K_p * pid_error[i] / D + flow_rate += factor * K_p * pid_error[i] / D end if !iszero(K_i) - flow_rate += reduction_factor * K_i * integral_value[i] / D + flow_rate += factor * K_i * integral_value[i] / D end if !iszero(K_d) @@ -751,11 +773,8 @@ function formulate_flow!( if active[i] hasindex, basin_idx = id_index(basin.node_id, upstream_basin_id) @assert hasindex "TabulatedRatingCurve must be downstream of a Basin" - s = storage[basin_idx] - reduction_factor = clamp(s, 0.0, 10.0) / 10.0 - q = - reduction_factor * - tables[i](get_level(p, upstream_basin_id, current_level, t)) + factor = reduction_factor(storage[basin_idx], 10.0) + q = factor * tables[i](get_level(p, upstream_basin_id, current_level, t)) else q = 0.0 end @@ -909,40 +928,84 @@ function formulate_flow!( end function formulate_flow!( - node::Union{Pump, Outlet}, + pump::Pump, p::Parameters, flow::AbstractMatrix, storage::AbstractVector, )::Nothing (; connectivity, basin) = p (; graph_flow) = connectivity - (; node_id, active, flow_rate, is_pid_controlled) = node + (; node_id, active, flow_rate, is_pid_controlled) = pump flow_rate = get_tmp(flow_rate, storage) for (id, isactive, rate, pid_controlled) in zip(node_id, active, flow_rate, is_pid_controlled) src_id = only(inneighbors(graph_flow, id)) dst_id = only(outneighbors(graph_flow, id)) - if !isactive + if !isactive || pid_controlled flow[src_id, id] = 0.0 flow[id, dst_id] = 0.0 continue end - if pid_controlled + hasindex, basin_idx = id_index(basin.node_id, src_id) + + q = rate + + if hasindex + # Pumping from basin + q *= reduction_factor(storage[basin_idx], 10.0) + end + + flow[src_id, id] = q + flow[id, dst_id] = q + end + return nothing +end + +function formulate_flow!( + outlet::Outlet, + p::Parameters, + flow::AbstractMatrix, + current_level::AbstractVector, + storage::AbstractVector, + t::Float64, +)::Nothing + (; connectivity, basin) = p + (; graph_flow) = connectivity + (; node_id, active, flow_rate, is_pid_controlled, min_crest_level) = outlet + flow_rate = get_tmp(flow_rate, storage) + for (i, id) in enumerate(node_id) + src_id = only(inneighbors(graph_flow, id)) + dst_id = only(outneighbors(graph_flow, id)) + + if !active[i] || is_pid_controlled[i] + flow[src_id, id] = 0.0 + flow[id, dst_id] = 0.0 continue end hasindex, basin_idx = id_index(basin.node_id, src_id) + q = flow_rate[i] + if hasindex - # Pumping from basin - s = storage[basin_idx] - reduction_factor = clamp(s, 0.0, 10.0) / 10.0 - q = reduction_factor * rate - else - # Pumping from level boundary - q = rate + # Flowing from basin + q *= reduction_factor(storage[basin_idx], 10.0) + end + + # No flow of outlet if source level is lower than target level + src_level = get_level(p, src_id, current_level, t) + dst_level = get_level(p, dst_id, current_level, t) + + if src_level !== nothing && dst_level !== nothing + Δlevel = src_level - dst_level + q *= reduction_factor(Δlevel, 0.1) + end + + # No flow out outlet if source level is lower than minimum crest level + if src_level !== nothing && !isnan(min_crest_level[i]) + q *= reduction_factor(src_level - min_crest_level[i], 0.1) end flow[src_id, id] = q @@ -995,7 +1058,7 @@ function formulate_flows!( formulate_flow!(flow_boundary, p, flow, t) formulate_flow!(fractional_flow, flow, p) formulate_flow!(pump, p, flow, storage) - formulate_flow!(outlet, p, flow, storage) + formulate_flow!(outlet, p, flow, current_level, storage, t) return nothing end @@ -1049,7 +1112,7 @@ function water_balance!( t, ) - # Negative storage musn't decrease, based on Shampine's et. al. advice + # Negative storage must not decrease, based on Shampine's et. al. advice # https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain for i in eachindex(u.storage) if u.storage[i] < 0 diff --git a/core/src/utils.jl b/core/src/utils.jl index d4b80baaa..65cb8372c 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -222,11 +222,11 @@ Ribasim.findlastgroup(2, [5,4,2,2,5,2,2,2,1]) """ function findlastgroup(id::Int, ids::AbstractVector{Int})::UnitRange{Int} idx_block_end = findlast(==(id), ids) - if isnothing(idx_block_end) + if idx_block_end === nothing return 1:0 end idx_block_begin = findprev(!=(id), ids, idx_block_end) - idx_block_begin = if isnothing(idx_block_begin) + idx_block_begin = if idx_block_begin === nothing 1 else # can happen if that id is the only ID in ids @@ -353,7 +353,7 @@ function set_static_value!( )::NamedTuple for (i, id) in enumerate(node_id) idx = findsorted(static.node_id, id) - isnothing(idx) && continue + idx === nothing && continue row = static[idx] set_table_row!(table, row, i) end @@ -380,7 +380,7 @@ function set_current_value!( row -> row.node_id == id && !isnan(getproperty(row, symbol)), pre_table, ) - if !isnothing(idx) + if idx !== nothing vector[i] = getproperty(pre_table, symbol)[idx] end end @@ -411,14 +411,18 @@ function get_level( node_id::Int, current_level::AbstractVector, t::Float64, -)::Real +)::Union{Real, Nothing} (; basin, level_boundary) = p hasindex, i = id_index(basin.node_id, node_id) return if hasindex current_level[i] else - i = findsorted(level_boundary.node_id, node_id)::Int - level_boundary.level[i](t) + i = findsorted(level_boundary.node_id, node_id) + if i === nothing + nothing + else + level_boundary.level[i](t) + end end end @@ -452,7 +456,7 @@ function basin_bottoms( )::Tuple{Float64, Float64} bottom_a = basin_bottom(basin, basin_a_id) bottom_b = basin_bottom(basin, basin_b_id) - if isnothing(bottom_a) && isnothing(bottom_b) + if bottom_a === bottom_b === nothing error(lazy"No bottom defined on either side of $id") end bottom_a = something(bottom_a, bottom_b) @@ -897,3 +901,18 @@ function Base.getindex(fv::FlatVector, i::Int) v = fv.v[d + 1] return v[r + 1] end + +""" +Function that goes smoothly from 0 to 1 in the interval [0,threshold], +and is constant outside this interval. +""" +function reduction_factor(x::T, threshold::Real)::T where {T <: Real} + return if x < 0 + zero(T) + elseif x < threshold + x_scaled = x / threshold + (-2 * x_scaled + 3) * x_scaled^2 + else + one(T) + end +end diff --git a/core/src/validation.jl b/core/src/validation.jl index 2eb1919d4..cafabd4db 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -141,6 +141,7 @@ end flow_rate::Float64 min_flow_rate::Union{Missing, Float64} max_flow_rate::Union{Missing, Float64} + min_crest_level::Union{Missing, Float64} control_state::Union{Missing, String} end diff --git a/core/test/control.jl b/core/test/control.jl index 4480099fe..51a30414a 100644 --- a/core/test/control.jl +++ b/core/test/control.jl @@ -188,6 +188,6 @@ end t_target_jump = discrete_control.record.time[2] t_idx_target_jump = searchsortedlast(timesteps, t_target_jump) - @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-4) - @test isapprox(level[end], target_low, atol = 1e-2) + @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-1) + @test isapprox(level[end], target_low, atol = 1e-1) end diff --git a/core/test/libribasim.jl b/core/test/libribasim.jl index 81ad77772..bb604b92d 100644 --- a/core/test/libribasim.jl +++ b/core/test/libribasim.jl @@ -19,9 +19,9 @@ toml_path = normpath(@__DIR__, "../../data/basic/basic.toml") toml_path_ptr = Base.unsafe_convert(Cstring, toml_path) # safe to finalize uninitialized model - @test isnothing(libribasim.model) + @test libribasim.model === nothing @test libribasim.finalize() == 0 - @test isnothing(libribasim.model) + @test libribasim.model === nothing # cannot get time of uninitialized model @test libribasim.last_error_message == "" @@ -41,6 +41,6 @@ toml_path = normpath(@__DIR__, "../../data/basic/basic.toml") @test unsafe_string(type_ptr) == "double" @test libribasim.finalize() == 0 - @test isnothing(libribasim.model) + @test libribasim.model === nothing end end diff --git a/core/test/run_models.jl b/core/test/run_models.jl index 1658824de..ad800deef 100644 --- a/core/test/run_models.jl +++ b/core/test/run_models.jl @@ -86,7 +86,7 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test successful_retcode(model) - @test model.integrator.sol.u[end] ≈ Float32[8.41143, 725.5068] skip = Sys.isapple() + @test model.integrator.sol.u[end] ≈ Float32[7.783636, 726.16394] skip = Sys.isapple() # the highest level in the dynamic table is updated to 1.2 from the callback @test model.integrator.p.tabulated_rating_curve.tables[end].t[end] == 1.2 end @@ -138,6 +138,33 @@ end @test A ≈ 10 * h end +@testset "Outlet constraints" begin + toml_path = normpath(@__DIR__, "../../data/outlet/outlet.toml") + @test ispath(toml_path) + + model = Ribasim.run(toml_path) + p = model.integrator.p + (; level_boundary, outlet) = p + (; level) = level_boundary + level = level[1] + + timesteps = model.saved_flow.t + outlet_flow = [saveval[1] for saveval in model.saved_flow.saveval] + + t_min_crest_level = + level.t[2] * (outlet.min_crest_level[1] - level.u[1]) / (level.u[2] - level.u[1]) + + # No outlet flow when upstream level is below minimum crest level + @test all(@. outlet_flow[timesteps <= t_min_crest_level] == 0) + + timesteps = Ribasim.timesteps(model) + t_maximum_level = level.t[2] + level_basin = Ribasim.get_storages_and_levels(model).level[:] + + # Basin level converges to stable level boundary level + all(isapprox.(level_basin[timesteps .>= t_maximum_level], level.u[3], atol = 5e-2)) +end + @testset "ManningResistance" begin """ Apply the "standard step method" finite difference method to find a diff --git a/core/test/utils.jl b/core/test/utils.jl index a0ae96bd7..0544cbe99 100644 --- a/core/test/utils.jl +++ b/core/test/utils.jl @@ -206,3 +206,13 @@ end @test isempty(fv) @test length(fv) == 0 end + +@testset "reduction_factor" begin + @test Ribasim.reduction_factor(-2.0, 2.0) === 0.0 + @test Ribasim.reduction_factor(0.0f0, 2.0) === 0.0f0 + @test Ribasim.reduction_factor(0.0, 2.0) === 0.0 + @test Ribasim.reduction_factor(1.0f0, 2.0) === 0.5f0 + @test Ribasim.reduction_factor(1.0, 2.0) === 0.5 + @test Ribasim.reduction_factor(3.0f0, 2.0) === 1.0f0 + @test Ribasim.reduction_factor(3.0, 2.0) === 1.0 +end diff --git a/core/test/validation.jl b/core/test/validation.jl index 64720ea41..f34c66184 100644 --- a/core/test/validation.jl +++ b/core/test/validation.jl @@ -230,6 +230,7 @@ end [-1.0], [NaN], [NaN], + [NaN], Dict{Tuple{Int, String}, NamedTuple}(), [false], ) diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index 1c29c6ffc..f2c884959 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -90,6 +90,62 @@ The presence of division by the basin area means that areas of size zero are not ::: # Natural water balance terms + +## The reduction factor {#sec-reduction_factor} +At several points in the equations below a *reduction factor* is used. This is a term that makes certain transitions more smooth, for instance when a pump stops providing water when its source basin dries up. The reduction factor is given by + +$$ + \phi(x; p) = + \begin{align} + \begin{cases} + 0 &\text{if}\quad x < 0 \\ + -2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\ + 1 &\text{if}\quad x > p + \end{cases} + \end{align}, +$$ + +where $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below. + +```{python} +# | code-fold: true +import numpy as np +import matplotlib.pyplot as plt + +def f(x, p = 3): + x_scaled = x / p + phi = (-2 * x_scaled + 3) * x_scaled**2 + phi = np.where(x < 0, 0, phi) + phi = np.where(x > p, 1, phi) + + return phi + +fontsize = 15 +p = 3 +N = 100 +x_min = -1 +x_max = 4 +x = np.linspace(x_min,x_max,N) +phi = f(x,p) + +fig,ax = plt.subplots(dpi=80) +ax.plot(x,phi) + +y_lim = ax.get_ylim() + +ax.set_xticks([0,p], [0,"$p$"], fontsize=fontsize) +ax.set_yticks([0,1], [0,1], fontsize=fontsize) +ax.hlines([0,1],x_min,x_max, color = "k", ls = ":", zorder=-1) +ax.vlines([0,p], *y_lim, color = "k", ls = ":") +ax.set_xlim(x_min,x_max) +ax.set_xlabel("$x$", fontsize=fontsize) +ax.set_ylabel("$\phi(x;p)$", fontsize=fontsize) +ax.set_ylim(y_lim) + +fig.tight_layout() +plt.show() +``` + ## Precipitation The precipitation term is given by @@ -107,14 +163,10 @@ surface may change as well, depending on the slopes of the surface waters. The evaporation term is given by $$ - Q_E = E_\text{pot} \cdot A(S) \cdot r. + Q_E = E_\text{pot} \cdot A(S) \cdot \phi(d;0.1). $$ {#eq-evap} -Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $r$ is the reduction factor which depends on the depth $d$ and is given by -$$ - r = \frac{\min\{d,0.1\}}{0.1} \le 1. -$$ {#eq-reduc} -It provides a smooth gradient as $S \rightarrow 0$ (but not at $d=0.1$). +Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $\phi$ is the [reduction factor](equations.qmd#sec-reduction_factor) which depends on the depth $d$. It provides a smooth gradient as $S \rightarrow 0$. A straightforward formulation $Q_E = \mathrm{max}(E_\text{pot} A(S), 0)$ is unsuitable, as $\frac{\mathrm{d}Q_E}{\mathrm{d}S}(S=0)$ is then not well-defined. @@ -183,26 +235,24 @@ basins (external nodes) in a network. 2. `LevelBoundary` 3. `FlowBoundary` -### Pump and outlet +### Pump {#sec-pump} -The behaviour of pumps and outlets is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate $q$, multiplied by a reduction factor: +The behaviour of pumps is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate $q$, multiplied by a reduction factor: $$ -Q_\text{pump/outlet} = \phi(u)q +Q_\text{pump} = \phi(u; 10.0)q $$ -Here $u$ is the storage of the upstream basin. The reduction factor +Here $u$ is the storage of the upstream basin. The [reduction factor](equations.qmd#sec-reduction_factor) $\phi$ makes sure that the flow of the pump goes smootly to $0$ as the upstream basin dries out. + +### Outlet {#sec-outlet} +The outlet is very similar to the pump, but it has a few extra [reduction factors](equations.qmd#sec-reduction_factor) for physical constraints: $$ -\phi(u) = -\begin{align} - \begin{cases} - \frac{u}{10} &\text{if}& 0 \leq u \le 10 \\ - 1 &\text{if}& u \geq 10 - \end{cases} -\end{align} -$$ {#eq-pumpreductionfactor} +Q_\text{outlet} = \phi(u_a; 10.0)\phi(\Delta h; 0.1) \phi(h_a-h_\text{min};0.1)q. +$$ +The subscript $a$ denotes the upstream node and $b$ the downstream node. The first reduction factor is equivalent to the one for the pump. The second one makes sure that the outlet flow goes to zero as the head difference $\Delta h = h_a - h_b$ goes to zero. The last one makes sure that the outlet only produces flow when the upstream level is above the minimum chrest level $h_\text{min}$. -makes sure that the flow of the pump goes to $0$ as the upstream basin dries out. +Not all node types upstream or downstream of the outlet have a defined level. If this is the case, and therefore the reduction factor cannot be computed, it is defined to be $1.0$. ### TabulatedRatingCurve @@ -415,19 +465,19 @@ $$ {#eq-PIDflow} for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so $$ -Q_\text{pump/outlet} = \text{clip}(\phi(u_\text{us})Q_\text{PID}; Q_\min, Q_\max). +Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_\min, Q_\max). $$ -Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\phi$ is the reduction factor (see @eq-pumpreductionfactor) and +Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\Phi$ is the product of [reduction factors](equations.qmd#sec-reduction_factor) associated with the [pump](equations.qmd#sec-pump) or [outlet](equations.qmd#sec-outlet) and $$ \text{clip}(Q; Q_\min, Q_\max) = \begin{align} \begin{cases} - Q_\min & \text{if} & Q < Q_\min, \\ - Q & \text{if} & Q_\min \leq Q \leq Q_\max, \\ - Q_\max & \text{if} & Q > Q_\max. - \end{cases} + Q_\min & \text{if} & Q < Q_\min \\ + Q & \text{if} & Q_\min \leq Q \leq Q_\max \\ + Q_\max & \text{if} & Q > Q_\max + \end{cases}. \end{align} $$ @@ -465,81 +515,6 @@ Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d $$ where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin's storage are known. -## Jacobian contributions - -For convenience here we denote a time derivative of a variable as a dot over the symbol. - -For the integral state we simply get the contribution -$$ -\frac{\partial \dot{I}}{\partial u_\text{PID}} = - \frac{1}{A(u_\text{PID})}. -$$ - -Note that when the calculated pump flow lies outside $[Q_\min, Q_\max]$, the pump flow is locally constant and thus all partial derivatives of $Q_\text{pump/outlet}$ with respect to states are $0$. Assuming that it lies in $[Q_\min, Q_\max]$, we derive the Jacobian contributions. Define the enumerator and denominator of the fraction in $Q_\text{pump/outlet}$: - -$$ -\begin{align} - E &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\ - D &= 1 \pm \phi(u_\text{us})\frac{K_d}{A(u_\text{PID})} -\end{align} -$$ - -Then - -$$ -\frac{\partial Q_\text{pump/outlet}}{\partial I} = K_i\frac{\phi(u_\text{us})}{D}. -$$ - -For the derivative of $Q_\text{pump/outlet}$ with respect to storages we need the derivatives of the enumerator and denominator. - -For the denominator we distinguish 2 possibilities: - -- If the controlled node is a pump, then $u_\text{us} \equiv u_\text{PID}$ and so $D$ is only a function of $u_\text{PID}$: -$$ -\frac{\text{d}D}{\text{d}u_\text{PID}} = - K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] -$$ -- If the controlled node is an outlet, then $D$ is a function of two distinct storages (assuming the upstream node is indeed a basin): -$$ -\begin{align} - \frac{\partial D}{\partial u_\text{PID}} &= \phi(u_\text{us})K_d\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}, \\ - \frac{\partial D}{\partial u_\text{us}} &= -\phi'(u_\text{us})K_d\frac{1}{A(u_\text{PID})}. -\end{align} -$$ - -For the enumerator there distinguish 2 cases: - -- The partial derivative with respect to $u_\text{PID}$: -$$ -\frac{\partial E}{\partial u_\text{PID}} = -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}. -$$ -- The partial derivative with respect to a different storage $u_i$: -$$ -\frac{\partial E}{\partial u_i} = - \frac{K_d}{A(u_\text{PID})}\frac{\partial \hat{f}_\text{PID}}{\partial u_i}. -$$ - -For computational efficiency we exploit that $\phi'$ is mostly $0$. Note that - -$$ -\frac{\partial\hat{f}_\text{PID}}{\partial u_i} = \hat{J}[i,\text{PID}], -$$ -\emph{i.e.} the Jacobian term of the dependence of $u_\text{PID}$ on $u_i$ without the contribution of the PID controlled pump/outlet. Note thus that also the Jacobian contribution of the PID controlled pump can only be computed after the other contributions to the Jacobian are known. - -The partial derivatives of $Q_\text{pump}$ are then given by -$$ -\begin{align} - \frac{\partial Q_\text{pump}}{\partial u_\text{PID}} &= \phi'(u_\text{PID})\frac{E}{D} + \frac{\phi(u_\text{PID})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\text{d}D}{\text{d}u_\text{PID}}\right], \\ - \frac{\partial Q_\text{pump}}{\partial u_i} &= \frac{\phi(u_\text{PID})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID}). -\end{align} -$$ - -The partial derivatives of $Q_\text{outlet}$ are given by -$$ -\begin{align} - \frac{\partial Q_\text{outlet}}{\partial u_\text{PID}} &= \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\partial D}{\partial u_\text{PID}}\right], \\ - \frac{\partial Q_\text{outlet}}{\partial u_\text{us}} &= \phi'(u_\text{us})\frac{E}{D} + \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{us}} - E\frac{\partial D}{\partial u_\text{us}}\right], \\ - \frac{\partial Q_\text{outlet}}{\partial u_i} &= \frac{\phi(u_\text{us})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID, us}). -\end{align} -$$ - ## The sign of the parameters Note by @eq-error that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level. diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 83afbdc67..526103d24 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -151,6 +151,8 @@ name it must have in the GeoPackage if it is stored there. - `TabulatedRatingCurve / time`: dynamic rating curve - Pump: pump water from a source node to a destination node - `Pump / static`: flow rate +- Outlet: let water flow with a prescribed flux under the conditions of positive head difference and the upstream level being higher than the minimum crest level + - `Outlet / static`: flow rate, minimum crest level - Terminal: Water sink without state or properties - `Terminal / static`: - (only node IDs) - DiscreteControl: Set parameters of other nodes based on model state conditions (e.g. basin level) @@ -309,12 +311,13 @@ level | Float64 | $m$ | - discharge | Float64 | $m^3 s^{-1}$ | non-negative -### Pump and outlet +### Pump -Pump/conduct water from a source node to a destination node. +Pump water from a source node to a destination node. The set flow rate will be pumped unless the intake storage is less than $10~m^3$, in which case the flow rate will be linearly reduced to $0~m^3/s$. The intake must be either a Basin or LevelBoundary. +When PID controlled, the pump must point away from the controlled basin in terms of edges. column | type | unit | restriction --------- | ------- | ------------ | ----------- @@ -325,6 +328,21 @@ min_flow_rate | Float64 | $m^3 s^{-1}$ | (optional, default 0.0) max_flow_rate | Float64 | $m^3 s^{-1}$ | (optional) control_state | String | - | (optional) +### Outlet + +The outlet is very similar to the pump. The outlet has two additional physical constraints: water only flows trough the outlet when the head difference is positive (i.e. water flows down by gravity), and the upstream level must be above the minimum crest level if the upstream level is defined. +When PID controlled, the outlet must point towards the controlled basin in terms of edges. + +column | type | unit | restriction +--------- | ------- | ------------ | ----------- +node_id | Int | - | sorted +active | Bool | - | (optional, default true) +flow_rate | Float64 | $m^3 s^{-1}$ | non-negative +min_flow_rate | Float64 | $m^3 s^{-1}$ | (optional, default 0.0) +max_flow_rate | Float64 | $m^3 s^{-1}$ | (optional) +min_crest_level | Float64 | $m$ | (optional) +control_state | String | - | (optional) + ### LevelBoundary Acts like an infinitely large basin where the level does not change by flow. diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index addae4e7a..acb00ca7b 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -635,7 +635,7 @@ "source": [ "# Model with discrete control\n", "\n", - "The model constructed below consists of a single basin which slowly drains trough a `TabulatedRatingCurve`, but is held within a range around a target level (setpoint) by two connected pumps." + "The model constructed below consists of a single basin which slowly drains trough a `TabulatedRatingCurve`, but is held within a range around a target level (setpoint) by two connected pumps. These two pumps behave like a reversible pump. When pumping can be done in only one direction, and the other direction is only possible under gravity, use an Outlet for that direction." ] }, { @@ -660,7 +660,7 @@ " (2.0, 0.0), # 4: LevelBoundary\n", " (-1.0, 0.0), # 5: TabulatedRatingCurve\n", " (-2.0, 0.0), # 6: Terminal\n", - " (0.0, 1.5), # 7: DiscreteControl\n", + " (1.0, 0.0), # 7: DiscreteControl\n", " ]\n", ")\n", "\n", @@ -1449,7 +1449,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/schema/OutletStatic.schema.json b/docs/schema/OutletStatic.schema.json index d307169fc..af03e5539 100644 --- a/docs/schema/OutletStatic.schema.json +++ b/docs/schema/OutletStatic.schema.json @@ -21,6 +21,13 @@ "boolean" ] }, + "min_crest_level": { + "format": "default", + "description": "min_crest_level", + "type": [ + "number" + ] + }, "flow_rate": { "format": "double", "description": "flow_rate", diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 998fdaeda..a61f9aa4d 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -147,7 +147,7 @@ def plot(self, ax=None, zorder=None) -> Any: "ManningResistance": "r", "TabulatedRatingCurve": "g", "Pump": "0.5", # grayscale level - "Outlet": "y", + "Outlet": "g", "Terminal": "m", "FlowBoundary": "m", "DiscreteControl": "k", diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 66c2e3bdc..4b271865f 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -200,6 +200,7 @@ def _write_tables(self, directory: FilePath) -> None: ) if isinstance(input_entry, TableModel) and not is_geometry: input_entry.write_table(connection) + connection.commit() shutil.move(temp_path, gpkg_path) return diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index da6f31259..9fb2ab244 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -1,6 +1,6 @@ # generated by datamodel-codegen: # filename: root.schema.json -# timestamp: 2023-08-21T14:05:19+00:00 +# timestamp: 2023-09-04T11:45:02+00:00 from __future__ import annotations @@ -130,6 +130,7 @@ class OutletStatic(BaseModel): max_flow_rate: Optional[float] = Field(None, description="max_flow_rate") remarks: Optional[str] = Field("", description="a hack for pandera") active: Optional[bool] = Field(None, description="active") + min_crest_level: Optional[float] = Field(None, description="min_crest_level") flow_rate: float = Field(..., description="flow_rate") node_id: int = Field(..., description="node_id") control_state: Optional[str] = Field(None, description="control_state") diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 49cad3bdf..2f15f4f77 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -4,6 +4,7 @@ from ribasim_testmodels.basic import ( basic_model, basic_transient_model, + outlet_model, tabulated_rating_curve_model, ) from ribasim_testmodels.bucket import bucket_model @@ -62,4 +63,5 @@ "invalid_edge_types_model", "discrete_control_of_pid_control_model", "level_boundary_condition_model", + "outlet_model", ] diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 8bbb72ccc..9da4892c2 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -5,7 +5,7 @@ def basic_model() -> ribasim.Model: - """Set up a basic model with all node types and static forcing""" + """Set up a basic model with most node types and static forcing""" # Setup the basins: profile = pd.DataFrame( @@ -379,3 +379,108 @@ def tabulated_rating_curve_model() -> ribasim.Model: ) return model + + +def outlet_model(): + """Set up a basic model with an outlet that encounters various physical constraints.""" + + # Set up the nodes: + xy = np.array( + [ + (0.0, 0.0), # 1: LevelBoundary + (1.0, 0.0), # 2: Outlet + (2.0, 0.0), # 3: Basin + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = ["LevelBoundary", "Outlet", "Basin"] + + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + static=gpd.GeoDataFrame( + data={"type": node_type}, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the edges: + from_id = np.array([1, 2], dtype=np.int64) + to_id = np.array([2, 3], dtype=np.int64) + lines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id) + edge = ribasim.Edge( + static=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": len(from_id) * ["flow"], + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": 3, + "area": 1000.0, + "level": [0.0, 1.0], + } + ) + + static = pd.DataFrame( + data={ + "node_id": [3], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + + state = pd.DataFrame(data={"node_id": [3], "level": 1e-3}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) + + # Setup the level boundary: + level_boundary = ribasim.LevelBoundary( + time=pd.DataFrame( + data={ + "node_id": 1, + "time": [ + "2020-01-01 00:00:00", + "2020-06-01 00:00:00", + "2021-01-01 00:00:00", + ], + "level": [1.0, 3.0, 3.0], + } + ) + ) + + # Setup the outlet + outlet = ribasim.Outlet( + static=pd.DataFrame( + data={ + "node_id": [2], + "flow_rate": 1e-3, + "min_crest_level": 2.0, + } + ) + ) + + model = ribasim.Model( + modelname="outlet", + node=node, + edge=edge, + basin=basin, + outlet=outlet, + level_boundary=level_boundary, + starttime="2020-01-01 00:00:00", + endtime="2021-01-01 00:00:00", + ) + + return model diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py index 7e7082435..b8fc718d9 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py @@ -113,7 +113,7 @@ def pid_control_model(): static=pd.DataFrame( data={ "node_id": [4], - "level": [1.0], # Not relevant + "level": [5.0], # Not relevant } ) ) @@ -290,7 +290,7 @@ def discrete_control_of_pid_control_model(): "node_id": [6, 6], "control_state": ["target_high", "target_low"], "listen_node_id": [3, 3], - "target": [6.0, 4.0], + "target": [5.0, 3.0], "proportional": 2 * [1e-2], "integral": 2 * [1e-8], "derivative": 2 * [-1e-1], diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py index 5f666e405..a12cb068a 100644 --- a/qgis/core/nodes.py +++ b/qgis/core/nodes.py @@ -177,7 +177,7 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "LevelBoundary": (QColor("green"), "LevelBoundary", shape.Circle), "FlowBoundary": (QColor("purple"), "FlowBoundary", shape.Hexagon), "Pump": (QColor("gray"), "Pump", shape.Hexagon), - "Outlet": (QColor("yellow"), "Outlet", shape.Hexagon), + "Outlet": (QColor("green"), "Outlet", shape.Hexagon), "ManningResistance": (QColor("red"), "ManningResistance", shape.Diamond), "Terminal": (QColor("purple"), "Terminal", shape.Square), "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), @@ -397,9 +397,10 @@ class LevelBoundaryStatic(Input): class LevelBoundaryTime(Input): - input_type = "LevelBoundary / static" + input_type = "LevelBoundary / time" geometry_type = "No Geometry" attributes = [ + QgsField("time", QVariant.DateTime), QgsField("node_id", QVariant.Int), QgsField("time", QVariant.DateTime), QgsField("level", QVariant.Double), @@ -428,6 +429,7 @@ class OutletStatic(Input): QgsField("flow_rate", QVariant.Double), QgsField("min_flow_rate", QVariant.Double), QgsField("max_flow_rate", QVariant.Double), + QgsField("min_crest_level", QVariant.Double), QgsField("control_state", QVariant.String), ]