Support transient user demand return factor

core/src/allocation_optim.jl

@@ -499,6 +499,7 @@ to the capacity of the outflow source.
 function adjust_capacities_returnflow!(
     allocation_model::AllocationModel,
     p::Parameters,
+    t::Float64,
 )::Nothing
     (; graph, user_demand) = p
     (; problem) = allocation_model
@@ -509,7 +510,7 @@ function adjust_capacities_returnflow!(
         constraint = constraints_outflow[node_id]
         capacity =
             JuMP.normalized_rhs(constraint) +
-            user_demand.return_factor[node_id.idx] *
+            user_demand.return_factor[node_id.idx](t) *
             JuMP.value(F[(inflow_id(graph, node_id), node_id)])
 
         JuMP.set_normalized_rhs(constraint, capacity)
@@ -973,7 +974,7 @@ function optimize_priority!(
     adjust_capacities_edge!(allocation_model)
     adjust_capacities_basin!(allocation_model)
     adjust_capacities_buffer!(allocation_model)
-    adjust_capacities_returnflow!(allocation_model, p)
+    adjust_capacities_returnflow!(allocation_model, p, t)
 
     # Adjust demands for next optimization (in case of internal_sources -> collect_demands)
     for parameter in propertynames(p)

core/src/parameter.jl

@@ -686,7 +686,7 @@ min_level: The level of the source basin below which the UserDemand does not abs
     demand_itp::Vector{Vector{ScalarInterpolation}}
     demand_from_timeseries::Vector{Bool}
     allocated::Matrix{Float64}
-    return_factor::Vector{Float64}
+    return_factor::Vector{ScalarInterpolation}
     min_level::Vector{Float64}
 end
 

core/src/read.jl

@@ -921,7 +921,7 @@
     active::Vector{Bool},
     demand::Matrix{Float64},
     demand_itp::Vector{Vector{ScalarInterpolation}},
-    return_factor::Vector{Float64},
+    return_factor::Vector{ScalarInterpolation},
     min_level::Vector{Float64},
     static::StructVector{UserDemandStaticV1},
     ids::Vector{Int32},
@@ -932,7 +932,12 @@
         user_demand_idx = searchsortedfirst(ids, first_row.node_id)
 
         active[user_demand_idx] = coalesce(first_row.active, true)
-        return_factor[user_demand_idx] = first_row.return_factor
+        return_factor_old = return_factor[user_demand_idx]
+        return_factor[user_demand_idx] = LinearInterpolation(
+            fill(first_row.return_factor, 2),
+            return_factor_old.t;
+            extrapolate = true,
+        )
         min_level[user_demand_idx] = first_row.min_level
 
         for row in group
@@ -955,7 +960,7 @@
     demand::Matrix{Float64},
     demand_itp::Vector{Vector{ScalarInterpolation}},
     demand_from_timeseries::Vector{Bool},
-    return_factor::Vector{Float64},
+    return_factor::Vector{ScalarInterpolation},
     min_level::Vector{Float64},
     time::StructVector{UserDemandTimeV1},
     ids::Vector{Int32},
@@ -971,11 +976,24 @@
 
         active[user_demand_idx] = true
         demand_from_timeseries[user_demand_idx] = true
-        return_factor[user_demand_idx] = first_row.return_factor
+        return_factor_itp, is_valid_return = get_scalar_interpolation(
+            config.starttime,
+            t_end,
+            StructVector(group),
+            NodeID(:UserDemand, first_row.node_id, 0),
+            :return_factor;
+        )
+        if is_valid_return
+            return_factor[user_demand_idx] = return_factor_itp
+        else
+            @error "The return_factor(t) relationship for UserDemand $(first_row.node_id) from the time table has repeated timestamps, this can not be interpolated."
+            errors = true
+        end
+
         min_level[user_demand_idx] = first_row.min_level
 
         priority_idx = findsorted(priorities, first_row.priority)
-        demand_p_itp, is_valid = get_scalar_interpolation(
+        demand_p_itp, is_valid_demand = get_scalar_interpolation(
             config.starttime,
             t_end,
             StructVector(group),
@@ -985,10 +1003,10 @@
         )
         demand[user_demand_idx, priority_idx] = demand_p_itp(0.0)
 
-        if is_valid
+        if is_valid_demand
             demand_itp[user_demand_idx][priority_idx] = demand_p_itp
         else
-            @error "The demand(t) relationship for UserDemand $node_id of priority $p from the time table has repeated timestamps, this can not be interpolated."
+            @error "The demand(t) relationship for UserDemand $(first_row.node_id) of priority $(first_row.priority_idx) from the time table has repeated timestamps, this can not be interpolated."
             errors = true
         end
     end
@@ -1022,7 +1040,8 @@
     ]
     demand_from_timeseries = fill(false, n_user)
     allocated = fill(Inf, n_user, n_priority)
-    return_factor = zeros(n_user)
+    return_factor =
+        [LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(node_ids)]
     min_level = zeros(n_user)
 
     # Process static table

core/src/solve.jl

@@ -273,7 +273,7 @@ function formulate_flow!(
         set_flow!(graph, inflow_edge, q)
 
         # Return flow is immediate
-        set_flow!(graph, outflow_edge, q * return_factor)
+        set_flow!(graph, outflow_edge, q * return_factor(t))
     end
     return nothing
 end

core/test/run_models_test.jl

@@ -385,18 +385,34 @@ end
 
 @testitem "UserDemand" begin
     using SciMLBase: successful_retcode
+    using Dates
+    using DataFrames: DataFrame
 
     toml_path = normpath(@__DIR__, "../../generated_testmodels/user_demand/ribasim.toml")
     @test ispath(toml_path)
     model = Ribasim.run(toml_path)
     @test successful_retcode(model)
 
-    day = 86400.0
-    @test only(model.integrator.sol(0day)) == 1000.0
-    # constant UserDemand withdraws to 0.9m/900m3
-    @test only(model.integrator.sol(150day)) ≈ 900 atol = 5
-    # dynamic UserDemand withdraws to 0.5m/509m3
-    @test only(model.integrator.sol(180day)) ≈ 509 atol = 1
+    seconds_in_day = 86400.0
+    @test only(model.integrator.sol(0seconds_in_day)) == 1000.0
+    # constant UserDemand withdraws to 0.9m or 900m3 due to min level = 0.9
+    @test only(model.integrator.sol(150seconds_in_day)) ≈ 900 atol = 5
+    # dynamic UserDemand withdraws to 0.5m or 500m3 due to min level = 0.5
+    @test only(model.integrator.sol(220seconds_in_day)) ≈ 500 atol = 1
+
+    # Trasient return factor
+    flow = DataFrame(Ribasim.flow_table(model))
+    return_factor_itp = model.integrator.p.user_demand.return_factor[3]
+    flow_in =
+        filter([:from_node_id, :to_node_id] => (from, to) -> (from, to) == (1, 4), flow)
+    flow_out =
+        filter([:from_node_id, :to_node_id] => (from, to) -> (from, to) == (4, 5), flow)
+    time_seconds = Ribasim.seconds_since.(flow_in.time, model.config.starttime)
+    @test isapprox(
+        flow_out.flow_rate,
+        return_factor_itp.(time_seconds) .* flow_in.flow_rate,
+        rtol = 1e-1,
+    )
 end
 
 @testitem "ManningResistance" begin

docs/guide/examples.ipynb

@@ -189,7 +189,7 @@
    "source": [
     "q = 10 / 86400  # 10 m³/day\n",
     "tabulated_rating_curve4 = model.tabulated_rating_curve.add(\n",
-    "    Node(8, Point(4.0, 0.0)),\n",
+    "    Node(8, Point(3.0, -1.0)),\n",
     "    [\n",
     "        tabulated_rating_curve.Static(\n",
     "            level=[0.0, 1.0],\n",
@@ -207,7 +207,7 @@
     "    ],\n",
     ")\n",
     "tabulated_rating_curve8 = model.tabulated_rating_curve.add(\n",
-    "    Node(4, Point(3.0, -1.0)),\n",
+    "    Node(4, Point(4.0, 0.0)),\n",
     "    [\n",
     "        tabulated_rating_curve.Static(\n",
     "            level=[0.0, 1.0],\n",
@@ -2339,7 +2339,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.4"
+   "version": "3.12.5"
   }
  },
  "nbformat": 4,

python/ribasim_api/tests/test_bmi.py

@@ -161,15 +161,15 @@ def test_get_value_ptr_allocation(libribasim, user_demand, tmp_path):
 
     # Demand
     actual_demand = libribasim.get_value_ptr("user_demand.demand")
-    expected_demand = np.array([1e-4, 0.0])
+    expected_demand = np.array([1e-4, 0.0, 0.0])
     assert_array_almost_equal(actual_demand, expected_demand)
 
     # Run model for a while to build up realized
     libribasim.update_until(60.0)
 
     # Realized
     actual_realized = libribasim.get_value_ptr("user_demand.realized")
-    expected_realized = np.array([1e-4 * 60.0, 0.0])
+    expected_realized = np.array([1e-4 * 60.0, 0.0, 0.0])
     assert_array_almost_equal(actual_realized, expected_realized)
 
 

python/ribasim_testmodels/ribasim_testmodels/allocation.py

@@ -59,12 +59,31 @@ def user_demand_model() -> Model:
             )
         ],
     )
-    model.terminal.add(Node(4, Point(2, 0)))
+    model.user_demand.add(
+        Node(4, Point(1, 0)),
+        [
+            user_demand.Time(
+                time=[
+                    "2020-08-01 00:00:00",
+                    "2020-09-01 00:00:00",
+                    "2020-10-01 00:00:00",
+                    "2020-11-01 00:00:00",
+                ],
+                min_level=0.0,
+                demand=[0.0, 1e-4, 2e-4, 0.0],
+                return_factor=[0.0, 0.1, 0.2, 0.3],
+                priority=1,
+            )
+        ],
+    )
+    model.terminal.add(Node(5, Point(2, 0)))
 
     model.edge.add(model.basin[1], model.user_demand[2])
     model.edge.add(model.basin[1], model.user_demand[3])
-    model.edge.add(model.user_demand[2], model.terminal[4])
-    model.edge.add(model.user_demand[3], model.terminal[4])
+    model.edge.add(model.basin[1], model.user_demand[4])
+    model.edge.add(model.user_demand[2], model.terminal[5])
+    model.edge.add(model.user_demand[3], model.terminal[5])
+    model.edge.add(model.user_demand[4], model.terminal[5])
 
     return model