Listen to concentration

core/src/callback.jl

@@ -314,7 +314,7 @@ Get a value for a condition. Currently supports getting levels from basins and f
 from flow boundaries.
 """
 function get_value(subvariable::NamedTuple, p::Parameters, u::AbstractVector, t::Float64)
-    (; flow_boundary, level_boundary) = p
+    (; flow_boundary, level_boundary, basin) = p
     (; listen_node_id, look_ahead, variable, variable_ref) = subvariable
 
     if !iszero(variable_ref.idx)
@@ -338,6 +338,8 @@ function get_value(subvariable::NamedTuple, p::Parameters, u::AbstractVector, t:
             error("Flow condition node $listen_node_id is not a flow boundary.")
         end
 
+    elseif startswith(variable, "concentration_external.")
+        value = basin.concentration_external[listen_node_id.idx][variable](t)
     else
         error("Unsupported condition variable $variable.")
     end

core/src/parameter.jl

@@ -268,6 +268,9 @@ end
     demand::Vector{Float64}
     # Data source for parameter updates
     time::StructVector{BasinTimeV1, C, Int}
+    # Concentrations
+    concentration_external::Vector{Dict{String, ScalarInterpolation}} =
+        Dict{String, ScalarInterpolation}[]
 end
 
 """

core/src/read.jl

@@ -584,6 +584,45 @@ function Basin(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int
     level_to_area = LinearInterpolation.(level_to_area)
     storage_to_level = invert_integral.(level_to_area)
 
+    t_end = seconds_since(config.endtime, config.starttime)
+
+    errors = false
+
+    concentration_external_data =
+        load_structvector(db, config, BasinConcentrationExternalV1)
+    concentration_external = Dict{String, ScalarInterpolation}[]
+    for id in node_id
+        concentration_external_id = Dict{String, ScalarInterpolation}()
+        data_id = filter(row -> row.node_id == id.value, concentration_external_data)
+        for group in IterTools.groupby(row -> row.substance, data_id)
+            first_row = first(group)
+            substance = first_row.substance
+            itp, no_duplication = get_scalar_interpolation(
+                config.starttime,
+                t_end,
+                StructVector(group),
+                NodeID(:Basin, first_row.node_id, 0),
+                :concentration,
+            )
+            concentration_external_id["concentration_external.$substance"] = itp
+            if any(itp.u .< 0)
+                errors = true
+                @error "Found negative concentration(s) in `Basin / concentration_external`." node_id =
+                    id, substance
+            end
+            if !no_duplication
+                errors = true
+                @error "There are repeated time values for in `Basin / concentration_external`." node_id =
+                    id substance
+            end
+        end
+        push!(concentration_external, concentration_external_id)
+    end
+
+    if errors
+        error("Errors encountered when parsing Basin concentration data.")
+    end
+
     return Basin(;
         node_id,
         inflow_ids = [collect(inflow_ids(graph, id)) for id in node_id],
@@ -599,6 +638,7 @@ function Basin(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int
         level_to_area,
         demand,
         time,
+        concentration_external,
     )
 end
 

core/src/schema.jl

@@ -54,7 +54,12 @@ function nodetype(
     record = Legolas.record_type(sv)
     node = last(split(string(Symbol(record)), '.'; limit = 3))
 
-    elements = split(string(T), '.'; limit = 3)
+    type_string = string(T)
+    elements = split(type_string, '.'; limit = 3)
+    last_element = last(elements)
+    if startswith(last_element, "concentration") && length(last_element) > 13
+        elements[end] = "concentration_$(last_element[14:end])"
+    end
     if isnode(sv)
         n = elements[2]
         k = Symbol(elements[3])

core/test/control_test.jl

@@ -281,3 +281,24 @@ end
     @test get_edge_flow("Basin", 1, "Outlet", 2) ≈ max.(0.4 .* inflow, 0) rtol = 1e-4
     @test get_edge_flow("Outlet", 2, "Terminal", 2) ≈ max.(0.4 .* inflow, 0) rtol = 1e-4
 end
+
+@testitem "Concentration discrete control" begin
+    using DataFrames: DataFrame
+
+    toml_path = normpath(
+        @__DIR__,
+        "../../generated_testmodels/concentration_condition/ribasim.toml",
+    )
+    @test ispath(toml_path)
+    model = Ribasim.run(toml_path)
+    flow_data = DataFrame(Ribasim.flow_table(model))
+    flow_edge_0 = filter(:edge_id => id -> id == 0, flow_data)
+    t = Ribasim.seconds_since.(flow_edge_0.time, model.config.starttime)
+    itp =
+        model.integrator.p.basin.concentration_external[1]["concentration_external.kryptonite"]
+    concentration = itp.(t)
+    threshold = 0.5
+    above_threshold = concentration .> threshold
+    @test all(isapprox.(flow_edge_0.flow_rate[above_threshold], 1e-3, rtol = 1e-2))
+    @test all(isapprox.(flow_edge_0.flow_rate[.!above_threshold], 0.0, atol = 1e-5))
+end

python/ribasim/ribasim/input_base.py

@@ -36,6 +36,22 @@
 __all__ = ("TableModel",)
 
 delimiter = " / "
+node_names_snake_case = [
+    "basin",
+    "continuous_control",
+    "discrete_control",
+    "flow_boundary",
+    "flow_demand",
+    "level_boundary",
+    "level_demand",
+    "linear_resistance",
+    "manning_resistance",
+    "outlet",
+    "pid_control",
+    "pump",
+    "tabulated_rating_curve",
+    "user_demand",
+]
 
 gpd.options.io_engine = "pyogrio"
 
@@ -181,11 +197,19 @@ def tablename(cls) -> str:
         NodeSchema -> Schema
         TabularRatingCurveStaticSchema -> TabularRatingCurve / Static
         """
-        names: list[str] = re.sub("([A-Z]+)", r" \1", str(cls.tableschema())).split()
-        if len(names) > 2:
-            return f"{''.join(names[:-2])}{delimiter}{names[-2].lower()}"
-        else:
+        cls_string = str(cls.tableschema())
+        names: list[str] = re.sub("([A-Z]+)", r" \1", cls_string).split()[:-1]
+        names_lowered = [name.lower() for name in names]
+        if len(names) == 1:
             return names[0]
+        else:
+            for n in range(1, len(names_lowered) + 1):
+                node_name_snake_case = "_".join(names_lowered[:n])
+                if node_name_snake_case in node_names_snake_case:
+                    node_name = "".join(names[:n])
+                    table_name = "_".join(names_lowered[n:])
+                    return node_name + delimiter + table_name
+            raise ValueError(f"Found no known node name in {cls_string}")
 
     @model_validator(mode="before")
     @classmethod

python/ribasim/ribasim/nodes/basin.py

@@ -15,13 +15,13 @@
 )
 
 __all__ = [
-    "Static",
-    "Time",
-    "State",
-    "Profile",
-    "Subgrid",
     "Area",
     "Concentration",
+    "Profile",
+    "State",
+    "Static",
+    "Subgrid",
+    "Time",
 ]
 
 

python/ribasim/tests/test_input_base.py

@@ -1,7 +1,19 @@
+from ribasim import geometry, nodes
 from ribasim.input_base import TableModel
 from ribasim.schemas import BasinSubgridSchema
 
 
 def test_tablemodel_schema():
     schema = TableModel[BasinSubgridSchema].tableschema()
     assert schema == BasinSubgridSchema
+
+
+def test_tablename():
+    cls = nodes.tabulated_rating_curve.Static
+    assert cls.tablename() == "TabulatedRatingCurve / static"
+
+    cls = nodes.basin.ConcentrationExternal
+    assert cls.tablename() == "Basin / concentration_external"
+
+    cls = geometry.edge.EdgeTable
+    assert cls.tablename() == "Edge"

python/ribasim_testmodels/ribasim_testmodels/__init__.py

@@ -30,6 +30,7 @@
 from ribasim_testmodels.continuous_control import outlet_continuous_control_model
 from ribasim_testmodels.discrete_control import (
     compound_variable_condition_model,
+    concentration_condition_model,
     connector_node_flow_condition_model,
     flow_condition_model,
     level_boundary_condition_model,
@@ -67,6 +68,7 @@
     "basic_transient_model",
     "bucket_model",
     "compound_variable_condition_model",
+    "concentration_condition_model",
     "connector_node_flow_condition_model",
     "discrete_control_of_pid_control_model",
     "fair_distribution_model",

python/ribasim_testmodels/ribasim_testmodels/discrete_control.py

@@ -1,3 +1,5 @@
+import numpy as np
+import pandas as pd
 from ribasim.config import Node
 from ribasim.model import Model
 from ribasim.nodes import (
@@ -528,3 +530,60 @@ def connector_node_flow_condition_model() -> Model:
     model.edge.add(model.discrete_control[1], model.linear_resistance[1])
 
     return model
+
+
+def concentration_condition_model() -> Model:
+    """
+    Set up a minimal model with discrete control based on a
+    concentration condition.
+    """
+
+    model = Model(
+        starttime="2020-01-01",
+        endtime="2021-01-01",
+        crs="EPSG:28992",
+    )
+
+    model.basin.add(
+        Node(1, Point(0, 0)),
+        [
+            basin.Profile(area=1000.0, level=[0.0, 1.0]),
+            basin.State(level=[20.0]),
+            basin.ConcentrationExternal(
+                time=pd.date_range(start="2020-01-01", end="2021-01-01", periods=100),
+                substance="kryptonite",
+                concentration=np.sin(np.linspace(0, 6 * np.pi, 100)) ** 2,
+            ),
+        ],
+    )
+
+    model.pump.add(
+        Node(1, Point(1, 0)),
+        [
+            pump.Static(
+                control_state=["On", "Off"], active=[True, False], flow_rate=1e-3
+            )
+        ],
+    )
+
+    model.terminal.add(Node(1, Point(2, 0)))
+
+    model.discrete_control.add(
+        Node(1, Point(1, 1)),
+        [
+            discrete_control.Variable(
+                listen_node_type=["Basin"],
+                listen_node_id=[1],
+                variable=["concentration_external.kryptonite"],
+                compound_variable_id=1,
+            ),
+            discrete_control.Condition(greater_than=[0.5], compound_variable_id=1),
+            discrete_control.Logic(truth_state=["T", "F"], control_state=["On", "Off"]),
+        ],
+    )
+
+    model.edge.add(model.basin[1], model.pump[1])
+    model.edge.add(model.pump[1], model.terminal[1])
+    model.edge.add(model.discrete_control[1], model.pump[1])
+
+    return model