Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support transient user demand return factor #1727

Merged
merged 10 commits into from
Aug 19, 2024
Merged

Support transient user demand return factor #1727

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
334bc32
use interpolation for return factor
Jingru923 Aug 14, 2024
281a963
address comments
Jingru923 Aug 14, 2024
54f5778
update a test model so transient return flow can be tested
Jingru923 Aug 14, 2024
24d05bc
update the unit tests
Jingru923 Aug 14, 2024
a4486e9
add one more user demand to bmi tests
Jingru923 Aug 14, 2024
98df891
new unit tests for return factor
Jingru923 Aug 15, 2024
e3bc943
Improve testing
SouthEndMusic Aug 18, 2024
cdd65b8
Unrelated: fix example plot
SouthEndMusic Aug 18, 2024
fb28c3f
Merge branch 'main' into return-factor
SouthEndMusic Aug 19, 2024
e24d2bf
final touch on the error message
Jingru923 Aug 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions core/src/allocation_optim.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion core/src/parameter.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
35 changes: 27 additions & 8 deletions core/src/read.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand All @@ -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
Expand All @@ -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},
Expand All @@ -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

Check warning on line 990 in core/src/read.jl

View check run for this annotation

Codecov / codecov/patch

core/src/read.jl#L989-L990 

Added lines #L989 - L990 were not covered by tests
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),
Expand All @@ -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 $priority_idx from the time table has repeated timestamps, this can not be interpolated."

Check warning on line 1009 in core/src/read.jl

View check run for this annotation

Codecov / codecov/patch

core/src/read.jl#L1009 

Added line #L1009 was not covered by tests
Jingru923 marked this conversation as resolved.
Show resolved Hide resolved
errors = true
end
end
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion core/src/solve.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions core/test/run_models_test.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -395,8 +395,8 @@ end
@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
# dynamic UserDemand withdraws to 0.5m/503m3
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is still too coarse a test. You should have a look at DataFrame(Ribasim.flow_table(model)), and check whether the outflow to inflow proportion matches your transient return factor

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Timepoint at 179day show a return ratio of 0.197ish and return factor is 0.2

@test only(model.integrator.sol(180day)) ≈ 503 atol = 1
end

@testitem "ManningResistance" begin
Expand Down
25 changes: 22 additions & 3 deletions python/ribasim_testmodels/ribasim_testmodels/allocation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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-06-01 00:00:00",
"2020-06-01 01:00:00",
"2020-07-01 00:00:00",
"2020-07-01 01:00:00",
],
demand=[0.0, 3e-4, 3e-4, 0.0],
return_factor=[0.0, 0.2, 0.3, 0.4],
min_level=0.5,
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

Expand Down