From f2473a27dbf71af8c99aacd98d058cf4e71a30aa Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 5 Sep 2023 17:02:34 +0200 Subject: [PATCH] `User` node type draft --- core/src/create.jl | 23 +++++++ core/src/solve.jl | 69 +++++++++++++++++++ core/src/validation.jl | 35 +++++++++- docs/contribute/addnode.qmd | 19 +++-- docs/schema/UserStatic.schema.json | 54 +++++++++++++++ docs/schema/UserTime.schema.json | 53 ++++++++++++++ docs/schema/root.schema.json | 6 ++ python/ribasim/ribasim/__init__.py | 2 + python/ribasim/ribasim/geometry/node.py | 2 + python/ribasim/ribasim/model.py | 4 ++ python/ribasim/ribasim/models.py | 24 ++++++- python/ribasim/ribasim/node_types/__init__.py | 2 + python/ribasim/ribasim/node_types/user.py | 46 +++++++++++++ qgis/core/nodes.py | 25 +++++++ 14 files changed, 353 insertions(+), 11 deletions(-) create mode 100644 docs/schema/UserStatic.schema.json create mode 100644 docs/schema/UserTime.schema.json create mode 100644 python/ribasim/ribasim/node_types/user.py diff --git a/core/src/create.jl b/core/src/create.jl index 590f7b600..6334db22f 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -621,6 +621,27 @@ function PidControl(db::DB, config::Config, chunk_size::Int)::PidControl ) end +function User(db::DB, config::Config)::User + static = load_structvector(db, config, UserStaticV1) + time = load_structvector(db, config, UserTimeV1) + parsed_parameters, valid = parse_static_and_time(db, config, "User"; static, time) + + 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.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 +659,7 @@ function Parameters(db::DB, config::Config)::Parameters 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 +678,7 @@ function Parameters(db::DB, config::Config)::Parameters terminal, discrete_control, pid_control, + user, Dict{Int, Symbol}(), ) for (fieldname, fieldtype) in zip(fieldnames(Parameters), fieldtypes(Parameters)) diff --git a/core/src/solve.jl b/core/src/solve.jl index 634acdab0..47a958222 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -403,6 +403,22 @@ struct PidControl{T} <: AbstractParameterNode control_mapping::Dict{Tuple{Int, String}, NamedTuple} end +""" +demand: water flux demand of user over time +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} + 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 +435,7 @@ struct Parameters{T, TSparse, C1, C2} terminal::Terminal discrete_control::DiscreteControl pid_control::PidControl{T} + user::User lookup::Dict{Int, Symbol} end @@ -719,6 +736,58 @@ function continuous_control!( return nothing end +function formulate!( + user::User, + p::Parameters, + current_level, + storage, + flow, + 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] + + # TODO: change to smooth reduction factor + # Smoothly let abstraction go to 0 as the source basin + # dries out + _, basin_idx = id_index(basin.node_id, src_id) + source_storage = storage[basin_idx] + reduction_factor_basin_empty = min(source_storage, 10.0) / 10.0 + q *= reduction_factor_basin_empty + + # TODO: change to smooth reduction factor + # 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] + reduction_factor_min_level = min(Δsource_level, 0.1) / 0.1 + q *= reduction_factor_min_level + + flow[src_id, id] = q + + # Return flow is immediate + flow[id, dst_id] = q * return_factor[i] + + return nothing + end +end + """ Directed graph: outflow is positive! """ diff --git a/core/src/validation.jl b/core/src/validation.jl index 4f422a5d6..505685ead 100644 --- a/core/src/validation.jl +++ b/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 @@ end 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{:Outlet}) = n_neighbor_bounds(1, 1, 1, typemax(Int) 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,7 +120,8 @@ n_neighbor_bounds_control(::Val{:Terminal}) = n_neighbor_bounds(0, 0, 0, 0) 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(nodetype) = +n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0) +n_neghbor_bounds_control(nodetype) = error("'n_neighbor_bounds_control' not defined for $nodetype.") # TODO NodeV1 and EdgeV1 are not yet used @@ -277,6 +288,24 @@ end control_state::Union{Missing, String} end +@version UserStaticV1 begin + node_id::Int + active::Union{Missing, String} + 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 diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd index a0924558c..b3ad1693f 100644 --- a/docs/contribute/addnode.qmd +++ b/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 diff --git a/docs/schema/UserStatic.schema.json b/docs/schema/UserStatic.schema.json new file mode 100644 index 000000000..4733934ca --- /dev/null +++ b/docs/schema/UserStatic.schema.json @@ -0,0 +1,54 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "properties": { + "remarks": { + "format": "default", + "default": "", + "description": "a hack for pandera", + "type": "string" + }, + "priority": { + "format": "default", + "description": "priority", + "type": "integer" + }, + "active": { + "format": "default", + "description": "active", + "type": [ + "string" + ] + }, + "demand": { + "format": "double", + "description": "demand", + "type": "number" + }, + "return_factor": { + "format": "default", + "description": "return_factor", + "type": "list" + }, + "node_id": { + "format": "default", + "description": "node_id", + "type": "integer" + }, + "allocated": { + "format": "double", + "description": "allocated", + "type": "number" + } + }, + "required": [ + "node_id", + "demand", + "allocated", + "return_factor", + "priority" + ], + "$id": "https://deltares.github.io/Ribasim/schema/UserStatic.schema.json", + "title": "UserStatic", + "description": "A UserStatic object based on Ribasim.UserStaticV1", + "type": "object" +} diff --git a/docs/schema/UserTime.schema.json b/docs/schema/UserTime.schema.json new file mode 100644 index 000000000..aa9eb10be --- /dev/null +++ b/docs/schema/UserTime.schema.json @@ -0,0 +1,53 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "properties": { + "remarks": { + "format": "default", + "default": "", + "description": "a hack for pandera", + "type": "string" + }, + "priority": { + "format": "default", + "description": "priority", + "type": "integer" + }, + "time": { + "format": "date-time", + "description": "time", + "type": "string" + }, + "demand": { + "format": "double", + "description": "demand", + "type": "number" + }, + "return_factor": { + "format": "default", + "description": "return_factor", + "type": "list" + }, + "node_id": { + "format": "default", + "description": "node_id", + "type": "integer" + }, + "allocated": { + "format": "double", + "description": "allocated", + "type": "number" + } + }, + "required": [ + "node_id", + "time", + "demand", + "allocated", + "return_factor", + "priority" + ], + "$id": "https://deltares.github.io/Ribasim/schema/UserTime.schema.json", + "title": "UserTime", + "description": "A UserTime object based on Ribasim.UserTimeV1", + "type": "object" +} diff --git a/docs/schema/root.schema.json b/docs/schema/root.schema.json index db2018158..469dcc354 100644 --- a/docs/schema/root.schema.json +++ b/docs/schema/root.schema.json @@ -10,12 +10,18 @@ "FlowBoundaryTime": { "$ref": "FlowBoundaryTime.schema.json" }, + "UserStatic": { + "$ref": "UserStatic.schema.json" + }, "PumpStatic": { "$ref": "PumpStatic.schema.json" }, "LevelBoundaryStatic": { "$ref": "LevelBoundaryStatic.schema.json" }, + "UserTime": { + "$ref": "UserTime.schema.json" + }, "DiscreteControlCondition": { "$ref": "DiscreteControlCondition.schema.json" }, diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index 7bf297f14..b02abbb88 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -18,6 +18,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.user import User __all__ = [ "models", @@ -40,4 +41,5 @@ "Terminal", "DiscreteControl", "PidControl", + "User", ] diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index a61f9aa4d..c604b781d 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -136,6 +136,7 @@ def plot(self, ax=None, zorder=None) -> Any: "FlowBoundary": "h", "DiscreteControl": "*", "PidControl": "x", + "User": "s", "": "o", } @@ -152,6 +153,7 @@ def plot(self, ax=None, zorder=None) -> Any: "FlowBoundary": "m", "DiscreteControl": "k", "PidControl": "k", + "User": "g", "": "k", } diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index aba0d84af..618017e05 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -33,6 +33,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.user import User from ribasim.types import FilePath @@ -75,6 +76,8 @@ class Model(BaseModel): Discrete control logic. pid_control : Optional[PidControl] PID controller attempting to set the level of a basin to a desired value using a pump/outlet. + user : Optional[User] + User node type with demand and priority. starttime : Union[str, datetime.datetime] Starting time of the simulation. endtime : Union[str, datetime.datetime] @@ -100,6 +103,7 @@ class Model(BaseModel): terminal: Optional[Terminal] discrete_control: Optional[DiscreteControl] pid_control: Optional[PidControl] + user: Optional[User] starttime: datetime.datetime endtime: datetime.datetime solver: Optional[Solver] diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index b9aabdfe2..e4ca81f28 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -4,7 +4,7 @@ from __future__ import annotations from datetime import datetime -from typing import Optional +from typing import Any, Optional from pydantic import BaseModel, Field @@ -31,6 +31,16 @@ class FlowBoundaryTime(BaseModel): node_id: int +class UserStatic(BaseModel): + remarks: Optional[str] = Field("", description="a hack for pandera") + priority: int = Field(..., description="priority") + active: Optional[str] = Field(None, description="active") + demand: float = Field(..., description="demand") + return_factor: Any = Field(..., description="return_factor") + node_id: int = Field(..., description="node_id") + allocated: float = Field(..., description="allocated") + + class PumpStatic(BaseModel): max_flow_rate: Optional[float] = None remarks: str = Field("", description="a hack for pandera") @@ -48,6 +58,16 @@ class LevelBoundaryStatic(BaseModel): level: float +class UserTime(BaseModel): + remarks: Optional[str] = Field("", description="a hack for pandera") + priority: int = Field(..., description="priority") + time: datetime = Field(..., description="time") + demand: float = Field(..., description="demand") + return_factor: Any = Field(..., description="return_factor") + node_id: int = Field(..., description="node_id") + allocated: float = Field(..., description="allocated") + + class DiscreteControlCondition(BaseModel): remarks: str = Field("", description="a hack for pandera") greater_than: float @@ -217,3 +237,5 @@ class Root(BaseModel): BasinProfile: Optional[BasinProfile] = None TerminalStatic: Optional[TerminalStatic] = None BasinStatic: Optional[BasinStatic] = None + UserStatic: Optional[UserStatic] = None + UserTime: Optional[UserTime] = None diff --git a/python/ribasim/ribasim/node_types/__init__.py b/python/ribasim/ribasim/node_types/__init__.py index f1671e7d5..a202a0883 100644 --- a/python/ribasim/ribasim/node_types/__init__.py +++ b/python/ribasim/ribasim/node_types/__init__.py @@ -10,6 +10,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.user import User __all__ = [ "Basin", @@ -24,4 +25,5 @@ "Terminal", "DiscreteControl", "PidControl", + "User", ] diff --git a/python/ribasim/ribasim/node_types/user.py b/python/ribasim/ribasim/node_types/user.py new file mode 100644 index 000000000..75e831d62 --- /dev/null +++ b/python/ribasim/ribasim/node_types/user.py @@ -0,0 +1,46 @@ +from typing import Optional + +import pandera as pa +from pandera.engines.pandas_engine import PydanticModel +from pandera.typing import DataFrame + +from ribasim import models +from ribasim.input_base import TableModel + +__all__ = ("User",) + + +class StaticSchema(pa.SchemaModel): + class Config: + """Config with dataframe-level data type.""" + + dtype = PydanticModel(models.UserStatic) + + +class TimeSchema(pa.SchemaModel): + class Config: + """Config with dataframe-level data type.""" + + dtype = PydanticModel(models.UserTime) + + +class User(TableModel): + """ + User node type with demand and priority. + + Parameters + ---------- + static: pandas.DataFrame + table with static data for this node type. + time: pandas.DataFrame + table with static data for this node type (only demand can be transient). + """ + + static: Optional[DataFrame[StaticSchema]] = None + time: Optional[DataFrame[TimeSchema]] = None + + class Config: + validate_assignment = True + + def sort(self): + self.static = self.static.sort_values("node_id", ignore_index=True) diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py index a12cb068a..fcb6dfdfc 100644 --- a/qgis/core/nodes.py +++ b/qgis/core/nodes.py @@ -182,6 +182,7 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "Terminal": (QColor("purple"), "Terminal", shape.Square), "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), "PidControl": (QColor("black"), "PidControl", shape.Cross2), + "User": (QColor("green"), "User", shape.Square), "": ( QColor("white"), "", @@ -509,6 +510,30 @@ class PidControlTime(Input): ] +class UserStatic(Input): + input_type = "User / static" + geometry_type = "No Geometry" + attributes = [ + QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), + QgsField("demand", QVariant.Double), + QgsField("return_factor", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + +class UserTime(Input): + input_type = "User / time" + geometry_type = "No Geometry" + attributes = [ + QgsField("node_id", QVariant.Int), + QgsField("time", QVariant.DateTime), + QgsField("demand", QVariant.Double), + QgsField("return_factor", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + NODES = {cls.input_type: cls for cls in Input.__subclasses__()} NONSPATIALNODETYPES = { cls.nodetype() for cls in Input.__subclasses__() if not cls.is_spatial()