User node type

core/src/bmi.jl

@@ -61,9 +61,13 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
             end
         end
 
-        # tstops for transient flow_boundary
+        # tell the solver to stop when new data comes in
+        # TODO add all time tables here
         time_flow_boundary = load_structvector(db, config, FlowBoundaryTimeV1)
-        tstops = get_tstops(time_flow_boundary.time, config.starttime)
+        tstops_flow_boundary = get_tstops(time_flow_boundary.time, config.starttime)
+        time_user = load_structvector(db, config, UserTimeV1)
+        tstops_user = get_tstops(time_user.time, config.starttime)
+        tstops = sort(unique(vcat(tstops_flow_boundary, tstops_user)))
 
         # use state
         state = load_structvector(db, config, BasinStateV1)

core/src/create.jl

@@ -411,7 +411,7 @@
 
 function Pump(db::DB, config::Config, chunk_size::Int)::Pump
     static = load_structvector(db, config, PumpStaticV1)
-    defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true)
+    defaults = (; min_flow_rate = 0.0, max_flow_rate = Inf, active = true)
     parsed_parameters, valid = parse_static_and_time(db, config, "Pump"; static, defaults)
     is_pid_controlled = falses(length(parsed_parameters.node_id))
 
@@ -440,7 +440,7 @@
 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, min_crest_level = NaN, active = true)
+        (; min_flow_rate = 0.0, max_flow_rate = Inf, min_crest_level = -Inf, active = true)
     parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults)
     is_pid_controlled = falses(length(parsed_parameters.node_id))
 
@@ -621,6 +621,38 @@
     )
 end
 
+function User(db::DB, config::Config)::User
+    static = load_structvector(db, config, UserStaticV1)
+    time = load_structvector(db, config, UserTimeV1)
+    defaults = (; min_level = -Inf, active = true)
+    time_interpolatables = [:demand]
+    parsed_parameters, valid = parse_static_and_time(
+        db,
+        config,
+        "User";
+        static,
+        time,
+        time_interpolatables,
+        defaults,
+    )
+
+    if !valid
+        error("Errors occurred when parsing User (node type) data.")
+    end
+
+    allocated = zeros(length(parsed_parameters.return_factor))
+
+    return User(
+        parsed_parameters.node_id,
+        parsed_parameters.active,
+        parsed_parameters.demand,
+        allocated,
+        parsed_parameters.return_factor,
+        parsed_parameters.min_level,
+        parsed_parameters.priority,
+    )
+end
+
 function Parameters(db::DB, config::Config)::Parameters
     n_states = length(get_ids(db, "Basin")) + length(get_ids(db, "PidControl"))
     chunk_size = pickchunksize(n_states)
@@ -638,6 +670,7 @@
     terminal = Terminal(db, config)
     discrete_control = DiscreteControl(db, config)
     pid_control = PidControl(db, config, chunk_size)
+    user = User(db, config)
 
     basin = Basin(db, config, chunk_size)
 
@@ -656,6 +689,7 @@
         terminal,
         discrete_control,
         pid_control,
+        user,
         Dict{Int, Symbol}(),
     )
     for (fieldname, fieldtype) in zip(fieldnames(Parameters), fieldtypes(Parameters))

core/src/solve.jl

@@ -403,6 +403,24 @@
     control_mapping::Dict{Tuple{Int, String}, NamedTuple}
 end
 
+"""
+demand: water flux demand of user over time
+active: whether this node is active and thus demands water
+allocated: water flux currently allocated to user
+return_factor: the factor in [0,1] of how much of the abstracted water is given back to the system
+min_level: The level of the source basin below which the user does not abstract
+priority: integer > 0, the lower the number the higher the priority of the users demand
+"""
+struct User
+    node_id::Vector{Int}
+    active::BitVector
+    demand::Vector{ScalarInterpolation}
+    allocated::Vector{Float64}
+    return_factor::Vector{Float64}
+    min_level::Vector{Float64}
+    priority::Vector{Int}
+end
+
 # TODO Automatically add all nodetypes here
 struct Parameters{T, TSparse, C1, C2}
     starttime::DateTime
@@ -419,6 +437,7 @@
     terminal::Terminal
     discrete_control::DiscreteControl
     pid_control::PidControl{T}
+    user::User
     lookup::Dict{Int, Symbol}
 end
 
@@ -719,6 +738,53 @@
     return nothing
 end
 
+function formulate_flow!(
+    user::User,
+    p::Parameters,
+    flow::AbstractMatrix,
+    current_level::AbstractVector,
+    storage::AbstractVector,
+    t::Float64,
+)::Nothing
+    (; connectivity, basin) = p
+    (; graph_flow) = connectivity
+    (; node_id, allocated, demand, active, return_factor, min_level) = user
+
+    for (i, id) in enumerate(node_id)
+        src_id = only(inneighbors(graph_flow, id))
+        dst_id = only(outneighbors(graph_flow, id))
+
+        if !active[i]
+            flow[src_id, id] = 0.0
+            flow[id, dst_id] = 0.0
+            continue
+        end
+
+        # For now allocated = demand
+        allocated[i] = demand[i](t)
+
+        q = allocated[i]
+
+        # Smoothly let abstraction go to 0 as the source basin dries out
+        _, basin_idx = id_index(basin.node_id, src_id)
+        factor_basin = reduction_factor(storage[basin_idx], 10.0)
+        q *= factor_basin
+
+        # Smoothly let abstraction go to 0 as the source basin
+        # level reaches its minimum level
+        source_level = get_level(p, src_id, current_level, t)
+        Δsource_level = source_level - min_level[i]
+        factor_level = reduction_factor(Δsource_level, 0.1)
+        q *= factor_level
+
+        flow[src_id, id] = q
+
+        # Return flow is immediate
+        flow[id, dst_id] = q * return_factor[i]
+    end
+    return nothing
+end
+
 """
 Directed graph: outflow is positive!
 """
@@ -1004,7 +1070,7 @@
         end
 
         # No flow out outlet if source level is lower than minimum crest level
-        if src_level !== nothing && !isnan(min_crest_level[i])
+        if src_level !== nothing
             q *= reduction_factor(src_level - min_crest_level[i], 0.1)
         end
 
@@ -1050,6 +1116,7 @@
         flow_boundary,
         pump,
         outlet,
+        user,
     ) = p
 
     formulate_flow!(linear_resistance, p, current_level, flow, t)
@@ -1059,6 +1126,7 @@
     formulate_flow!(fractional_flow, flow, p)
     formulate_flow!(pump, p, flow, storage)
     formulate_flow!(outlet, p, flow, current_level, storage, t)
+    formulate_flow!(user, p, flow, current_level, storage, t)
 
     return nothing
 end

core/src/validation.jl

@@ -22,6 +22,8 @@
 @schema "ribasim.tabulatedratingcurve.static" TabulatedRatingCurveStatic
 @schema "ribasim.tabulatedratingcurve.time" TabulatedRatingCurveTime
 @schema "ribasim.outlet.static" OutletStatic
+@schema "ribasim.user.static" UserStatic
+@schema "ribasim.user.time" UserTime
 
 const delimiter = " / "
 tablename(sv::Type{SchemaVersion{T, N}}) where {T, N} = join(nodetype(sv), delimiter)
@@ -46,8 +48,15 @@
 neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype))
 neighbortypes(::Val{:Pump}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary))
 neighbortypes(::Val{:Outlet}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary))
-neighbortypes(::Val{:Basin}) =
-    Set((:LinearResistance, :TabulatedRatingCurve, :ManningResistance, :Pump, :Outlet))
+neighbortypes(::Val{:User}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary))
+neighbortypes(::Val{:Basin}) = Set((
+    :LinearResistance,
+    :TabulatedRatingCurve,
+    :ManningResistance,
+    :Pump,
+    :Outlet,
+    :User,
+))
 neighbortypes(::Val{:Terminal}) = Set{Symbol}() # only endnode
 neighbortypes(::Val{:FractionalFlow}) = Set((:Basin, :Terminal, :LevelBoundary))
 neighbortypes(::Val{:FlowBoundary}) =
@@ -93,6 +102,7 @@
 n_neighbor_bounds_flow(::Val{:Terminal}) = n_neighbor_bounds(1, typemax(Int), 0, 0)
 n_neighbor_bounds_flow(::Val{:PidControl}) = n_neighbor_bounds(0, 0, 0, 0)
 n_neighbor_bounds_flow(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 0, 0)
+n_neighbor_bounds_flow(::Val{:User}) = n_neighbor_bounds(1, 1, 1, 1)
 n_neighbor_bounds_flow(nodetype) =
     error("'n_neighbor_bounds_flow' not defined for $nodetype.")
 
@@ -110,6 +120,7 @@
 n_neighbor_bounds_control(::Val{:PidControl}) = n_neighbor_bounds(0, 1, 1, 1)
 n_neighbor_bounds_control(::Val{:DiscreteControl}) =
     n_neighbor_bounds(0, 0, 1, typemax(Int))
+n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0)
 n_neighbor_bounds_control(nodetype) =
     error("'n_neighbor_bounds_control' not defined for $nodetype.")
 
@@ -277,6 +288,24 @@
     control_state::Union{Missing, String}
 end
 
+@version UserStaticV1 begin
+    node_id::Int
+    active::Union{Missing, Bool}
+    demand::Float64
+    return_factor::Float64
+    min_level::Float64
+    priority::Int
+end
+
+@version UserTimeV1 begin
+    node_id::Int
+    time::DateTime
+    demand::Float64
+    return_factor::Float64
+    min_level::Float64
+    priority::Int
+end
+
 function variable_names(s::Any)
     filter(x -> !(x in (:node_id, :control_state)), fieldnames(s))
 end
@@ -327,6 +356,7 @@
     LevelBoundaryTimeV1,
     PidControlTimeV1,
     TabulatedRatingCurveTimeV1,
+    UserTimeV1,
 }
 
 function sort_by_function(table::StructVector{<:TimeSchemas})

core/test/run_models.jl

@@ -165,6 +165,19 @@ end
     all(isapprox.(level_basin[timesteps .>= t_maximum_level], level.u[3], atol = 5e-2))
 end
 
+@testset "User" begin
+    toml_path = normpath(@__DIR__, "../../data/user/user.toml")
+    @test ispath(toml_path)
+    model = Ribasim.run(toml_path)
+
+    day = 86400.0
+    @test only(model.integrator.sol(0day)) == 1000.0
+    # constant user withdraws to 0.9m/900m3
+    @test only(model.integrator.sol(150day)) ≈ 900 atol = 5
+    # dynamic user withdraws to 0.5m/500m3
+    @test only(model.integrator.sol(180day)) ≈ 500 atol = 1
+end
+
 @testset "ManningResistance" begin
     """
     Apply the "standard step method" finite difference method to find a

docs/contribute/addnode.qmd

@@ -47,15 +47,19 @@ Now we define the function that is called in the second bullet above, in `create
 ```julia
 function NewNodeType(db::DB, config::Config)::NewNodeType
     static = load_structvector(db, config, NewNodeTypeStaticV1)
-    defaults = (; active = true)
+    defaults = (; foo = 1, bar = false)
     # Process potential control states in the static data
-    static_parsed = parse_static(static, db, "Outlet", defaults)
+    parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults)
+
+    if !valid
+        error("Errors occurred when parsing NewNodeType data.")
+    end
 
     # Unpack the fields of static as inputs for the NewNodeType constructor
     return NewNodeType(
-        static_parsed.node_id,
-        static_parsed.some_property,
-        static_parsed.control_mapping)
+        parsed_parameters.node_id,
+        parsed_parameters.some_property,
+        parsed_parameters.control_mapping)
 end
 ```
 
@@ -93,7 +97,7 @@ The current dependency groups are:
 - Either-neighbor dependencies: examples are `LinearResistance`, `ManningResistance`. If either the in-neighbor or out-neighbor of a node of this group is a basin, the storage of this basin depends on itself. If both the in-neighbor and the out-neighbor are basins, their storages also depend on eachother.
 - The `PidControl` node is a special case which is discussed in [equations](../core/equations.qmd#sec-PID).
 
-In the methods `formulate_jac!` in `jac.jl` the analytical expressions for the partial derivatives $\frac{\partial Q_{i',j'}}{\partial u_i}$ (and the ones related to PID integral states) are hardcoded. For `NewNodeType` either a new method of `formulate_jac!` has to be introduced, or it has to be added to the list of node types that do not contribute to the Jacobian in the method of `formulate_jac!` whose signature contains `node::AbstractParameterNode`.
+Using `jac_prototype` the Jacobian of `water_balance!` is computed automatically using [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/) with memory management provided by [PreallocationTools.jl](https://docs.sciml.ai/PreallocationTools/stable/). These computations make use of `DiffCache` and dual numbers.
 
 # Python I/O
 
@@ -102,7 +106,8 @@ In the methods `formulate_jac!` in `jac.jl` the analytical expressions for the p
 Create a new file `python/ribasim/ribasim/node_types/new_node_type.py` which is structured as follows:
 
 ```python
-import pandas as pd
+from typing import Optional
+
 import pandera as pa
 from pandera.engines.pandas_engine import PydanticModel
 from pandera.typing import DataFrame
@@ -139,10 +144,8 @@ class NewNodeType(TableModel):
     class Config:
         validate_assignment = True
 
-
-
     def sort(self):
-        self.static = self.static.sort_values("node_id", ignore_index=True)
+        self.static.sort_values("node_id", ignore_index=True, inplace=True)
 
 ```
 The `sort` method should implement the same sorting as in `validation.jl`.