diff --git a/Manifest.toml b/Manifest.toml index d03ba3a69..425ca14c5 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "d694251e41cfe037d14e0d129b2e2e49035c0990" +project_hash = "f8b2dea16bf0ccef5eef2111258ff045598662c9" [[deps.ADTypes]] git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0" @@ -375,9 +375,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[deps.DiffEqBase]] deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"] -git-tree-sha1 = "d7d43a1cc11dc4e4e5378816ae720fccd75a77c8" +git-tree-sha1 = "fa7d580038451a8df4434ef5b079ac9b2d486194" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.155.0" +version = "6.155.1" [deps.DiffEqBase.extensions] DiffEqBaseCUDAExt = "CUDA" @@ -1344,7 +1344,7 @@ uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" version = "3.5.18" [[deps.Ribasim]] -deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqCallbacks", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] +deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqBase", "DiffEqCallbacks", "EnumX", "FiniteDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] path = "core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" version = "2024.10.0" diff --git a/Project.toml b/Project.toml index f8fd65409..ef7bcf793 100644 --- a/Project.toml +++ b/Project.toml @@ -15,11 +15,11 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" diff --git a/core/Project.toml b/core/Project.toml index 521d86e91..76be2b5cd 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -15,10 +15,10 @@ DBInterface = "a10d1c49-ce27-4219-8d33-6db1a4562965" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" @@ -67,10 +67,10 @@ DataFrames = "1.4" DataInterpolations = "6" DataStructures = "0.18" Dates = "<0.0.1, 1" +DiffEqBase = "6.155" DiffEqCallbacks = "3.6" EnumX = "1.0" FiniteDiff = "2.21" -ForwardDiff = "0.10" Graphs = "1.9" HiGHS = "1.7" IOCapture = "0.2" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 06c1bf7c3..0c0ab18f5 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -21,6 +21,7 @@ using OrdinaryDiffEqCore: get_du, AbstractNLSolver, calculate_residuals! +using DiffEqBase: DiffEqBase using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs! using LineSearches: BackTracking diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 915a0139e..dc9613688 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -306,12 +306,12 @@ function get_basin_data( u::ComponentVector, node_id::NodeID, ) - (; graph, allocation) = p + (; graph, allocation, basin) = p (; Δt_allocation) = allocation_model (; mean_input_flows) = allocation @assert node_id.type == NodeType.Basin influx = mean_input_flows[(node_id, node_id)][] - storage_basin = u.storage[node_id.idx] + storage_basin = basin.current_storage[parent(u)][node_id.idx] control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(control_inneighbors) level_demand_idx = 0 diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 0c11a5dbe..0c12c251e 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -32,7 +32,7 @@ end function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{Float64} if name == "basin.storage" - model.integrator.u.storage + model.integrator.p.basin.current_storage[parent(model.integrator.u)] elseif name == "basin.level" model.integrator.p.basin.current_level[parent(model.integrator.u)] elseif name == "basin.infiltration" diff --git a/core/src/callback.jl b/core/src/callback.jl index be5a08d7b..f0a4393c3 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -7,21 +7,32 @@ Returns the CallbackSet and the SavedValues for flow. function create_callbacks( parameters::Parameters, config::Config, + u0::ComponentVector, saveat, )::Tuple{CallbackSet, SavedResults} (; starttime, basin, tabulated_rating_curve) = parameters callbacks = SciMLBase.DECallback[] + # Check for negative storage negative_storage_cb = FunctionCallingCallback(check_negative_storage) push!(callbacks, negative_storage_cb) - integrating_flows_cb = FunctionCallingCallback(integrate_flows!; func_start = false) - push!(callbacks, integrating_flows_cb) + # Save storages and levels + saved_basin_states = SavedValues(Float64, SavedBasinState) + save_basin_state_cb = SavingCallback(save_basin_state, saved_basin_states; saveat) + push!(callbacks, save_basin_state_cb) + # Update cumulative flows (exact integration and for BMI) + integrated_flows_cb = + FunctionCallingCallback(update_cumulative_flows!; func_start = false) + push!(callbacks, integrated_flows_cb) + + # Update Basin forcings tstops = get_tstops(basin.time.time, starttime) basin_cb = PresetTimeCallback(tstops, update_basin!; save_positions = (false, false)) push!(callbacks, basin_cb) + # Update TabulatedRatingCurve Q(h) relationships tstops = get_tstops(tabulated_rating_curve.time.time, starttime) tabulated_rating_curve_cb = PresetTimeCallback( tstops, @@ -33,13 +44,9 @@ function create_callbacks( # If saveat is a vector which contains 0.0 this callback will still be called # at t = 0.0 despite save_start = false saveat = saveat isa Vector ? filter(x -> x != 0.0, saveat) : saveat - saved_vertical_flux = SavedValues(Float64, typeof(basin.vertical_flux_integrated)) - save_vertical_flux_cb = - SavingCallback(save_vertical_flux, saved_vertical_flux; saveat, save_start = false) - push!(callbacks, save_vertical_flux_cb) - # save the flows over time - saved_flow = SavedValues(Float64, SavedFlow) + # save the flows averaged over the saveat intervals + saved_flow = SavedValues(Float64, SavedFlow{typeof(u0)}) save_flow_cb = SavingCallback(save_flow, saved_flow; saveat, save_start = false) push!(callbacks, save_flow_cb) @@ -66,7 +73,7 @@ function create_callbacks( saved = SavedResults( saved_flow, - saved_vertical_flux, + saved_basin_states, saved_subgrid_level, saved_solver_stats, ) @@ -75,152 +82,276 @@ function create_callbacks( return callback, saved end -function save_solver_stats(u, t, integrator) - (; stats) = integrator.sol - (; - time = t, - rhs_calls = stats.nf, - linear_solves = stats.nsolve, - accepted_timesteps = stats.naccept, - rejected_timesteps = stats.nreject, - ) -end - -function check_negative_storage(u, t, integrator)::Nothing - (; basin) = integrator.p - (; node_id) = basin - errors = false - for id in node_id - if u.storage[id.idx] < 0 - @error "Negative storage detected in $id" - errors = true - end - end - - if errors - t_datetime = datetime_since(integrator.t, integrator.p.starttime) - error("Negative storages found at $t_datetime.") - end - return nothing -end - """ -Integrate flows over the last timestep +Update with the latest timestep: +- Cumulative flows/forcings which are exposed via the BMI +- Cumulative flows/forcings which are integrated exactly +- Cumulative flows/forcings which are input for the allocation algorithm +- Cumulative flows/forcings which are realized demands in the allocation context + """ -function integrate_flows!(u, t, integrator)::Nothing - (; p, dt) = integrator - (; graph, user_demand, basin, allocation) = p - (; flow, flow_dict, flow_prev, flow_integrated) = graph[] - (; vertical_flux, vertical_flux_prev, vertical_flux_integrated, vertical_flux_bmi) = - basin - du = get_du(integrator) - water_balance!(du, u, p, t) - flow = flow[parent(du)] - vertical_flux = vertical_flux[parent(du)] - if !isempty(flow_prev) && isnan(flow_prev[1]) - # If flow_prev is not populated yet - copyto!(flow_prev, flow) +function update_cumulative_flows!(u, t, integrator)::Nothing + (; p, uprev, tprev, dt) = integrator + (; basin, user_demand, flow_boundary, allocation) = p + (; vertical_flux_bmi, vertical_flux_from_input) = basin + + # Update tprev + p.tprev[] = t + + # Update cumulative flows exposed via the BMI + @. user_demand.realized_bmi += u.user_demand_inflow - uprev.user_demand_inflow + @. vertical_flux_bmi.drainage += vertical_flux_from_input.drainage * dt + @. vertical_flux_bmi.evaporation += u.evaporation - uprev.evaporation + @. vertical_flux_bmi.infiltration += u.infiltration - uprev.infiltration + + # Update cumulative forcings which are integrated exactly + @. basin.cumulative_drainage += vertical_flux_from_input.drainage * dt + @. basin.cumulative_drainage_saveat += vertical_flux_from_input.drainage * dt + + # Precipitation depends on fixed area + for node_id in basin.node_id + fixed_area = basin_areas(basin, node_id.idx)[end] + added_precipitation = + fixed_area * vertical_flux_from_input.precipitation[node_id.idx] * dt + + vertical_flux_bmi.precipitation[node_id.idx] += added_precipitation + basin.cumulative_precipitation[node_id.idx] += added_precipitation + basin.cumulative_precipitation_saveat[node_id.idx] += added_precipitation end - @. flow_integrated += 0.5 * (flow + flow_prev) * dt - @. vertical_flux_integrated += 0.5 * (vertical_flux + vertical_flux_prev) * dt - @. vertical_flux_bmi += 0.5 * (vertical_flux + vertical_flux_prev) * dt + # Exact boundary flow over time step + for (id, flow_rate, active) in + zip(flow_boundary.node_id, flow_boundary.flow_rate, flow_boundary.active) + if active + volume = integral(flow_rate, tprev, t) + flow_boundary.cumulative_flow[id.idx] += volume + flow_boundary.cumulative_flow_saveat[id.idx] += volume + end + end - # UserDemand realized flows for BMI - for id in user_demand.node_id - src_id = inflow_id(graph, id) - flow_idx = flow_dict[src_id, id] - user_demand.realized_bmi[id.idx] += - 0.5 * (flow[flow_idx] + flow_prev[flow_idx]) * dt + # Update realized flows for allocation input + for edge in keys(allocation.mean_input_flows) + allocation.mean_input_flows[edge] += flow_update_on_edge(integrator, edge) end - # *Demand realized flow for output - for (edge, value) in allocation.mean_realized_flows - if edge[1] !== edge[2] - value += - 0.5 * - (get_flow(graph, edge..., du) + get_flow_prev(graph, edge..., du)) * - dt - allocation.mean_realized_flows[edge] = value + # Update realized flows for allocation output + for edge in keys(allocation.mean_realized_flows) + allocation.mean_realized_flows[edge] += flow_update_on_edge(integrator, edge) + if edge[1] == edge[2] + basin_id = edge[1] + @assert basin_id.type == NodeType.Basin + for inflow_id in basin.inflow_ids[basin_id.idx] + allocation.mean_realized_flows[edge] += + flow_update_on_edge(integrator, (inflow_id, basin_id)) + end + for outflow_id in basin.outflow_ids[basin_id.idx] + allocation.mean_realized_flows[edge] -= + flow_update_on_edge(integrator, (basin_id, outflow_id)) + end end end + return nothing +end - # Allocation source flows - for (edge, value) in allocation.mean_input_flows - if edge[1] == edge[2] - # Vertical fluxes - allocation.mean_input_flows[edge] = - value + - 0.5 * - ( - get_influx(basin, edge[1].idx) + - get_influx(basin, edge[1].idx; prev = true) - ) * - dt +""" +Given an edge (from_id, to_id), compute the cumulative flow over that +edge over the latest timestep. If from_id and to_id are both the same basin, +the function returns the sum of the basin forcings. +""" +function flow_update_on_edge( + integrator::DEIntegrator, + edge_src::Tuple{NodeID, NodeID}, +)::Float64 + (; u, uprev, p, t, tprev, dt) = integrator + (; basin, flow_boundary) = p + (; vertical_flux_from_input) = basin + from_id, to_id = edge_src + if from_id == to_id + @assert from_id.type == to_id.type == NodeType.Basin + idx = from_id.idx + fixed_area = basin_areas(basin, idx)[end] + ( + fixed_area * vertical_flux_from_input.precipitation[idx] + + vertical_flux_from_input.drainage[idx] + ) * dt - (u.evaporation[idx] - uprev.evaporation[idx]) - + (u.infiltration[idx] - uprev.infiltration[idx]) + elseif from_id.type == NodeType.FlowBoundary + if flow_boundary.active[from_id.idx] + integral(flow_boundary.flow_rate[from_id.idx], tprev, t) else - # Horizontal flows - allocation.mean_input_flows[edge] = - value + - 0.5 * - (get_flow(graph, edge..., du) + get_flow_prev(graph, edge..., du)) * - dt + 0.0 end + else + flow_idx = get_state_index(u, edge_src) + u[flow_idx] - uprev[flow_idx] end - copyto!(flow_prev, flow) - copyto!(vertical_flux_prev, vertical_flux) - return nothing end -"Compute the average flows over the last saveat interval and write -them to SavedValues" -function save_flow(u, t, integrator) - (; graph) = integrator.p - (; flow_integrated, flow_dict) = graph[] - (; node_id) = integrator.p.basin +""" +Save the storages and levels at the latest t. +""" +function save_basin_state(u, t, integrator) + (; p) = integrator + (; basin) = p + du = get_du(integrator) + current_storage = basin.current_storage[parent(du)] + current_level = basin.current_level[parent(du)] + water_balance!(du, u, p, t) + SavedBasinState(; storage = copy(current_storage), level = copy(current_level)) +end +""" +Save all cumulative forcings and flows over edges over the latest timestep, +Both computed by the solver and integrated exactly. Also computes the total horizontal +inflow and outflow per basin. +""" +function save_flow(u, t, integrator) + (; p, sol) = integrator + (; basin, state_inflow_edge, state_outflow_edge, flow_boundary) = p Δt = get_Δt(integrator) - flow_mean = copy(flow_integrated) - flow_mean ./= Δt - fill!(flow_integrated, 0.0) - - # Divide the flows over edges to Basin inflow and outflow, regardless of edge direction. - inflow_mean = zeros(length(node_id)) - outflow_mean = zeros(length(node_id)) - - for basin_id in node_id - for inflow_id in inflow_ids(graph, basin_id) - q = flow_mean[flow_dict[inflow_id, basin_id]] - if q > 0 - inflow_mean[basin_id.idx] += q + flow_mean = (u - sol(t - Δt)) / Δt + + inflow_mean = zeros(length(basin.node_id)) + outflow_mean = zeros(length(basin.node_id)) + + # Flow contributions from horizontal flow states + for (flow, inflow_edge, outflow_edge) in + zip(flow_mean, state_inflow_edge, state_outflow_edge) + inflow_id = inflow_edge.edge[1] + if inflow_id.type == NodeType.Basin + if flow > 0 + outflow_mean[inflow_id.idx] += flow else - outflow_mean[basin_id.idx] -= q + inflow_mean[inflow_id.idx] -= flow end end - for outflow_id in outflow_ids(graph, basin_id) - q = flow_mean[flow_dict[basin_id, outflow_id]] - if q > 0 - outflow_mean[basin_id.idx] += q + + outflow_id = outflow_edge.edge[2] + if outflow_id.type == NodeType.Basin + if flow > 0 + inflow_mean[outflow_id.idx] += flow else - inflow_mean[basin_id.idx] -= q + outflow_mean[outflow_id.idx] -= flow + end + end + end + + # Flow contributions from flow boundaries + flow_boundary_mean = copy(flow_boundary.cumulative_flow_saveat) ./ Δt + flow_boundary.cumulative_flow_saveat .= 0.0 + + for (outflow_edges, id) in zip(flow_boundary.outflow_edges, flow_boundary.node_id) + flow = flow_boundary_mean[id.idx] + for outflow_edge in outflow_edges + outflow_id = outflow_edge.edge[2] + if outflow_id.type == NodeType.Basin + inflow_mean[outflow_id.idx] += flow end end end - return SavedFlow(; flow = flow_mean, inflow = inflow_mean, outflow = outflow_mean) + precipitation = copy(basin.cumulative_precipitation_saveat) ./ Δt + drainage = copy(basin.cumulative_drainage_saveat) ./ Δt + @. basin.cumulative_precipitation_saveat = 0.0 + @. basin.cumulative_drainage_saveat = 0.0 + + saved_flow = SavedFlow(; + flow = flow_mean, + inflow = inflow_mean, + outflow = outflow_mean, + flow_boundary = flow_boundary_mean, + precipitation, + drainage, + ) + check_water_balance_error(integrator, saved_flow, Δt) + return saved_flow end -"Compute the average vertical fluxes over the last saveat interval and write -them to SavedValues" -function save_vertical_flux(u, t, integrator) +function check_water_balance_error( + integrator::DEIntegrator, + saved_flow::SavedFlow, + Δt::Float64, +)::Nothing + (; u, p, t) = integrator + (; basin, water_balance_abstol, water_balance_reltol) = p + errors = false + current_storage = basin.current_storage[parent(u)] + formulate_storages!(current_storage, u, u, p, t) + + for ( + inflow_rate, + outflow_rate, + precipitation, + drainage, + evaporation, + infiltration, + s_now, + s_prev, + id, + ) in zip( + saved_flow.inflow, + saved_flow.outflow, + saved_flow.precipitation, + saved_flow.drainage, + saved_flow.flow.evaporation, + saved_flow.flow.infiltration, + current_storage, + basin.storage_prev_saveat, + basin.node_id, + ) + storage_rate = (s_now - s_prev) / Δt + total_in = inflow_rate + precipitation + drainage + total_out = outflow_rate + evaporation + infiltration + balance_error = storage_rate - (total_in - total_out) + mean_flow_rate = (total_in + total_out) / 2 + relative_error = iszero(mean_flow_rate) ? 0.0 : balance_error / mean_flow_rate + + if abs(balance_error) > water_balance_abstol && + abs(relative_error) > water_balance_reltol + errors = true + @error "Too large water balance error" id balance_error relative_error + end + end + if errors + t = datetime_since(t, p.starttime) + error("Too large water balance error(s) detected at t = $t") + end + + @. basin.storage_prev_saveat = current_storage + return nothing +end + +function save_solver_stats(u, t, integrator) + (; stats) = integrator.sol + (; + time = t, + rhs_calls = stats.nf, + linear_solves = stats.nsolve, + accepted_timesteps = stats.naccept, + rejected_timesteps = stats.nreject, + ) +end + +function check_negative_storage(u, t, integrator)::Nothing (; basin) = integrator.p - (; vertical_flux_integrated) = basin + (; node_id, current_storage) = basin + du = get_du(integrator) + set_current_basin_properties!(du, u, integrator.p, t) + current_storage = current_storage[parent(du)] - Δt = get_Δt(integrator) - vertical_flux_mean = copy(vertical_flux_integrated) - vertical_flux_mean ./= Δt - fill!(vertical_flux_integrated, 0.0) + errors = false + for id in node_id + if current_storage[id.idx] < 0 + @error "Negative storage detected in $id" + errors = true + end + end - return vertical_flux_mean + if errors + t_datetime = datetime_since(integrator.t, integrator.p.starttime) + error("Negative storages found at $t_datetime.") + end + return nothing end """ @@ -243,6 +374,7 @@ function apply_discrete_control!(u, t, integrator)::Nothing (; discrete_control) = p (; node_id) = discrete_control du = get_du(integrator) + water_balance!(du, u, p, t) # Loop over the discrete control nodes to determine their truth state # and detect possible control state changes @@ -436,12 +568,10 @@ end "Load updates from 'Basin / time' into the parameters" function update_basin!(integrator)::Nothing - (; p, u) = integrator + (; p) = integrator (; basin) = p - (; storage) = u - (; node_id, time, vertical_flux_from_input, vertical_flux, vertical_flux_prev) = basin + (; node_id, time, vertical_flux_from_input) = basin t = datetime_since(integrator.t, integrator.p.starttime) - vertical_flux = vertical_flux[parent(u)] rows = searchsorted(time.time, t) timeblock = view(time, rows) @@ -457,41 +587,37 @@ function update_basin!(integrator)::Nothing i = searchsortedfirst(node_id, NodeID(NodeType.Basin, row.node_id, 0)) set_table_row!(table, row, i) end - - update_vertical_flux!(basin, storage) - - # Forget about vertical fluxes to handle discontinuous forcing from basin_update - copyto!(vertical_flux_prev, vertical_flux) return nothing end "Solve the allocation problem for all demands and assign allocated abstractions." function update_allocation!(integrator)::Nothing (; p, t, u) = integrator - (; allocation) = p + (; allocation, basin) = p + (; current_storage) = basin (; allocation_models, mean_input_flows, mean_realized_flows) = allocation + # Make sure current storages are up to date + du = get_du(integrator) + current_storage = current_storage[parent(du)] + formulate_storages!(current_storage, du, u, p, t) + # Don't run the allocation algorithm if allocation is not active # (Specifically for running Ribasim via the BMI) if !is_active(allocation) return nothing end - (; Δt_allocation) = allocation_models[1] - # Divide by the allocation Δt to obtain the mean input flows # from the integrated flows - for key in keys(mean_input_flows) - mean_input_flows[key] /= Δt_allocation + (; Δt_allocation) = allocation_models[1] + for edge in keys(mean_input_flows) + mean_input_flows[edge] /= Δt_allocation end # Divide by the allocation Δt to obtain the mean realized flows # from the integrated flows - for (edge, value) in mean_realized_flows - if edge[1] == edge[2] - # Compute the mean realized demand for basins as Δstorage/Δt_allocation - mean_realized_flows[edge] = value + u[edge[1].idx] - end + for edge in keys(mean_realized_flows) mean_realized_flows[edge] /= Δt_allocation end @@ -515,13 +641,6 @@ function update_allocation!(integrator)::Nothing mean_flows[edge] = 0.0 end end - - # Set basin storages for mean storage change computation - for (edge, value) in mean_realized_flows - if edge[1] == edge[2] - mean_realized_flows[edge] = value - u[edge[1].idx] - end - end end "Load updates from 'TabulatedRatingCurve / time' into the parameters" diff --git a/core/src/config.jl b/core/src/config.jl index 0d7b624d5..95e4d8c52 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -102,6 +102,8 @@ const nodetypes = collect(keys(nodekinds)) force_dtmin::Bool = false abstol::Float64 = 1e-6 reltol::Float64 = 1e-5 + water_balance_abstol::Float64 = 1e-3 + water_balance_reltol::Float64 = 1e-2 maxiters::Int = 1e9 sparse::Bool = true autodiff::Bool = true diff --git a/core/src/graph.jl b/core/src/graph.jl index f0efacf83..d948511ec 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -30,11 +30,10 @@ function create_graph(db::DB, config::Config)::MetaGraph node_ids = Dict{Int32, Set{NodeID}}() # Source edges per subnetwork edges_source = Dict{Int32, Set{EdgeMetadata}}() - # The flow counter gives a unique consecutive id to the - # flow edges to index the flow vectors - flow_counter = 0 + # The metadata of the flow edges in the order in which they are in the input + # and will be in the output + flow_edges = EdgeMetadata[] # Dictionary from flow edge to index in flow vector - flow_dict = Dict{Tuple{NodeID, NodeID}, Int}() graph = MetaGraph( DiGraph(); label_type = NodeID, @@ -80,20 +79,18 @@ function create_graph(db::DB, config::Config)::MetaGraph end edge_metadata = EdgeMetadata(; id = edge_id, - flow_idx = edge_type == EdgeType.flow ? flow_counter + 1 : 0, type = edge_type, subnetwork_id_source = subnetwork_id, edge = (id_src, id_dst), ) + if edge_type == EdgeType.flow + push!(flow_edges, edge_metadata) + end if haskey(graph, id_src, id_dst) errors = true @error "Duplicate edge" id_src id_dst end graph[id_src, id_dst] = edge_metadata - if edge_type == EdgeType.flow - flow_counter += 1 - flow_dict[(id_src, id_dst)] = flow_counter - end if subnetwork_id != 0 if !haskey(edges_source, subnetwork_id) edges_source[subnetwork_id] = Set{EdgeMetadata}() @@ -109,20 +106,7 @@ function create_graph(db::DB, config::Config)::MetaGraph error("Incomplete connectivity in subnetwork") end - flow = cache(flow_counter) - flow_prev = fill(NaN, flow_counter) - flow_integrated = zeros(flow_counter) - flow_edges = [edge for edge in values(graph.edge_data) if edge.type == EdgeType.flow] - graph_data = (; - node_ids, - edges_source, - flow_edges, - flow_dict, - flow, - flow_prev, - flow_integrated, - config.solver.saveat, - ) + graph_data = (; node_ids, edges_source, flow_edges, config.solver.saveat) graph = @set graph.graph_data = graph_data return graph @@ -183,48 +167,6 @@ function Base.iterate(iter::OutNeighbors, state = 1) return label_out, state end -function set_flow!(graph::MetaGraph, edge_metadata::EdgeMetadata, q::Number, du)::Nothing - set_flow!(graph, edge_metadata.flow_idx, q, du) - return nothing -end - -function set_flow!(graph, flow_idx::Int, q::Number, du)::Nothing - (; flow) = graph[] - flow[parent(du)][flow_idx] = q - return nothing -end - -""" -Get the flow over the given edge (du is needed for LazyBufferCache from ForwardDiff.jl). -""" -function get_flow(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, du)::Number - (; flow_dict) = graph[] - flow_idx = flow_dict[id_src, id_dst] - return get_flow(graph, flow_idx, du) -end - -function get_flow(graph, edge_metadata::EdgeMetadata, du)::Number - return get_flow(graph, edge_metadata.flow_idx, du) -end - -function get_flow(graph::MetaGraph, flow_idx::Int, du) - return graph[].flow[parent(du)][flow_idx] -end - -function get_flow_prev(graph, id_src::NodeID, id_dst::NodeID, du)::Number - (; flow_dict) = graph[] - flow_idx = flow_dict[id_src, id_dst] - return get_flow(graph, flow_idx, du) -end - -function get_flow_prev(graph, edge_metadata::EdgeMetadata, du)::Number - return get_flow_prev(graph, edge_metadata.flow_idx, du) -end - -function get_flow_prev(graph::MetaGraph, flow_idx::Int, du) - return graph[].flow_prev[parent(du)][flow_idx] -end - """ Get the inneighbor node IDs of the given node ID (label) over the given edge type in the graph. @@ -300,11 +242,38 @@ function inflow_id(graph::MetaGraph, id::NodeID)::NodeID end """ -Get the metadata of an edge in the graph from an edge of the underlying -DiGraph. +Get the specific q from the input vector `flow` which has the same components as +the state vector, given an edge (inflow_id, outflow_id). +`flow` can be either instantaneous or integrated/averaged. Instantaneous FlowBoundary flows can be obtained +from the parameters, but integrated/averaged FlowBoundary flows must be provided via `boundary_flow`. """ -function metadata_from_edge(graph::MetaGraph, edge::Edge{Int})::EdgeMetadata - label_src = label_for(graph, edge.src) - label_dst = label_for(graph, edge.dst) - return graph[label_src, label_dst] +function get_flow( + flow::ComponentVector, + p::Parameters, + t::Number, + edge::Tuple{NodeID, NodeID}; + boundary_flow = nothing, +) + (; flow_boundary) = p + from_id = edge[1] + if from_id.type == NodeType.FlowBoundary + if boundary_flow === nothing + flow_boundary.active[from_id.idx] ? flow_boundary.flow_rate[from_id.idx](t) : + 0.0 + else + boundary_flow[from_id.idx] + end + else + flow[get_state_index(flow, edge)] + end +end + +function get_influx(du::ComponentVector, id::NodeID, p::Parameters) + @assert id.type == NodeType.Basin + (; basin) = p + (; vertical_flux_from_input) = basin + fixed_area = basin_areas(basin, id.idx)[end] + return fixed_area * vertical_flux_from_input.precipitation[id.idx] + + vertical_flux_from_input.drainage[id.idx] - du.evaporation[id.idx] - + du.infiltration[id.idx] end diff --git a/core/src/logging.jl b/core/src/logging.jl index 4ee7617fa..517b227a8 100644 --- a/core/src/logging.jl +++ b/core/src/logging.jl @@ -35,12 +35,12 @@ function log_bottlenecks(model; converged::Bool) # Indicate convergence bottlenecks if possible with the current algorithm if hasproperty(cache, :nlsolver) - storage_error = @. abs(cache.nlsolver.cache.atmp.storage / u.storage) - perm = sortperm(storage_error; rev = true) + flow_error = @. abs(cache.nlsolver.cache.atmp / u) + perm = sortperm(flow_error; rev = true) errors = Pair{Symbol, Float64}[] for i in perm - node_id = Symbol(basin.node_id[i]) - error = storage_error[i] + node_id = Symbol(id_from_state_index(p, u, i)) + error = flow_error[i] if error < model.config.solver.reltol break end diff --git a/core/src/model.jl b/core/src/model.jl index 8b932c521..cec9e599d 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -1,6 +1,6 @@ -struct SavedResults{V1 <: ComponentVector{Float64}} - flow::SavedValues{Float64, SavedFlow} - vertical_flux::SavedValues{Float64, V1} +struct SavedResults{V <: ComponentVector{Float64}} + flow::SavedValues{Float64, SavedFlow{V}} + basin_state::SavedValues{Float64, SavedBasinState} subgrid_level::SavedValues{Float64, Vector{Float64}} solver_stats::SavedValues{Float64, SolverStats} end @@ -94,36 +94,29 @@ function Model(config::Config)::Model push!(tstops, get_tstops(time_schema.time, config.starttime)) end - # use state - state = load_structvector(db, config, BasinStateV1) - n = length(get_ids(db, "Basin")) - finally # always close the database, also in case of an error close(db) end @debug "Read database into memory." - storage = get_storages_from_levels(parameters.basin, state.level) - - @assert length(storage) == n "Basin / state length differs from number of Basins" - # Integrals for PID control - integral = zeros(length(parameters.pid_control.node_id)) - u0 = ComponentVector{Float64}(; storage, integral) + u0 = build_state_vector(parameters) du0 = zero(u0) + parameters = set_state_flow_edges(parameters, u0) + # The Solver algorithm alg = algorithm(config.solver; u0) - # Synchronize level with storage - set_current_basin_properties!(parameters.basin, u0, du0) - # for Float32 this method allows max ~1000 year simulations without accuracy issues t_end = seconds_since(config.endtime, config.starttime) @assert eps(t_end) < 3600 "Simulation time too long" t0 = zero(t_end) timespan = (t0, t_end) + # Synchronize level with storage + set_current_basin_properties!(du0, u0, parameters, t0) + saveat = convert_saveat(config.solver.saveat, t_end) saveat isa Float64 && push!(tstops, range(0, t_end; step = saveat)) tstops = sort(unique(vcat(tstops...))) @@ -139,7 +132,7 @@ function Model(config::Config)::Model prob = ODEProblem(RHS, u0, timespan, parameters) @debug "Setup ODEProblem." - callback, saved = create_callbacks(parameters, config, saveat) + callback, saved = create_callbacks(parameters, config, u0, saveat) @debug "Created callbacks." # Run water_balance! before initializing the integrator. This is because @@ -159,7 +152,7 @@ function Model(config::Config)::Model progress_steps = 100, callback, tstops, - isoutofdomain = (u, p, t) -> any(<(0), u.storage), + isoutofdomain, saveat, adaptive, dt, @@ -215,7 +208,6 @@ function SciMLBase.step!(model::Model, dt::Float64)::Model if round(ntimes) ≈ ntimes update_allocation!(integrator) end - set_previous_flows!(integrator) step!(integrator, dt, true) return model end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 2399003ab..73ca33f59 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -213,7 +213,6 @@ end """ Type for storing metadata of edges in the graph: id: ID of the edge (only used for labeling flow output) -flow_idx: Index in the vector of flows type: type of the edge subnetwork_id_source: ID of subnetwork where this edge is a source (0 if not a source) @@ -221,12 +220,13 @@ edge: (from node ID, to node ID) """ @kwdef struct EdgeMetadata id::Int32 - flow_idx::Int type::EdgeType.T subnetwork_id_source::Int32 edge::Tuple{NodeID, NodeID} end +Base.length(::EdgeMetadata) = 1 + """ The update of an parameter given by a value and a reference to the target location of the variable in memory @@ -258,14 +258,28 @@ abstract type AbstractDemandNode <: AbstractParameterNode end """ In-memory storage of saved mean flows for writing to results. -- `flow`: The mean flows on all edges +- `flow`: The mean flows on all edges and state-dependent forcings - `inflow`: The sum of the mean flows coming into each basin - `outflow`: The sum of the mean flows going out of each basin +- `flow_boundary`: The exact integrated mean flows of flow boundaries +- `precipitation`: The exact integrated mean precipitation +- `drainage`: The exact integrated mean drainage """ -@kwdef struct SavedFlow - flow::Vector{Float64} +@kwdef struct SavedFlow{V} + flow::V inflow::Vector{Float64} outflow::Vector{Float64} + flow_boundary::Vector{Float64} + precipitation::Vector{Float64} + drainage::Vector{Float64} +end + +""" +In-memory storage of saved instantaneous storages and levels for writing to results. +""" +@kwdef struct SavedBasinState + storage::Vector{Float64} + level::Vector{Float64} end """ @@ -290,13 +304,21 @@ end outflow_ids::Vector{Vector{NodeID}} = [NodeID[]] # Vertical fluxes vertical_flux_from_input::V1 = zeros(length(node_id)) - vertical_flux::Cache = cache(length(node_id)) - vertical_flux_prev::V2 = zeros(length(node_id)) - vertical_flux_integrated::V2 = zeros(length(node_id)) vertical_flux_bmi::V2 = zeros(length(node_id)) + # Initial_storage + storage0::Vector{Float64} = zeros(length(node_id)) + storage_prev_saveat::Vector{Float64} = zeros(length(node_id)) + # Analytically integrated forcings + cumulative_precipitation::Vector{Float64} = zeros(length(node_id)) + cumulative_drainage::Vector{Float64} = zeros(length(node_id)) + cumulative_precipitation_saveat::Vector{Float64} = zeros(length(node_id)) + cumulative_drainage_saveat::Vector{Float64} = zeros(length(node_id)) # Cache this to avoid recomputation + current_storage::Cache = cache(length(node_id)) current_level::Cache = cache(length(node_id)) current_area::Cache = cache(length(node_id)) + current_cumulative_precipitation::Cache = cache(length(node_id)) + current_cumulative_drainage::Cache = cache(length(node_id)) # Discrete values for interpolation storage_to_level::Vector{ LinearInterpolationIntInv{ @@ -437,12 +459,16 @@ end node_id: node ID of the FlowBoundary node outflow_edges: The outgoing flow edge metadata active: whether this node is active and thus contributes flow -flow_rate: target flow rate +cumulative_flow: The exactly integrated cumulative boundary flow since the start of the simulation +cumulative_flow_saveat: The exactly integrated cumulative boundary flow since the last saveat +flow_rate: flow rate (exact) """ @kwdef struct FlowBoundary <: AbstractParameterNode node_id::Vector{NodeID} outflow_edges::Vector{Vector{EdgeMetadata}} active::Vector{Bool} + cumulative_flow::Vector{Float64} = zeros(length(node_id)) + cumulative_flow_saveat::Vector{Float64} = zeros(length(node_id)) flow_rate::Vector{ScalarInterpolation} end @@ -578,15 +604,21 @@ node_id: node ID of the Terminal node end """ -A variant on `Base.Ref` where the source array is a vector that is possibly wrapped in a ForwardDiff.DiffCache. +A variant on `Base.Ref` where the source array is a vector that is possibly wrapped in a ForwardDiff.LazyBufferCache, +or a reference to the state derivative vector du. Retrieve value with get_value(ref::PreallocationRef, val) where `val` determines the return type. """ struct PreallocationRef vector::Cache idx::Int + from_du::Bool + function PreallocationRef(vector::Cache, idx::Int; from_du = false) + new(vector, idx, from_du) + end end -get_value(ref::PreallocationRef, du) = ref.vector[parent(du)][ref.idx] +get_value(ref::PreallocationRef, du) = + ref.from_du ? du[ref.idx] : ref.vector[parent(du)][ref.idx] function set_value!(ref::PreallocationRef, value, du)::Nothing ref.vector[parent(du)][ref.idx] = value @@ -760,12 +792,7 @@ node_ids: mapping subnetwork ID -> node IDs in that subnetwork edges_source: mapping subnetwork ID -> metadata of allocation source edges in that subnetwork flow_edges: The metadata of all flow edges -flow dict: mapping (source ID, destination ID) -> index in the flow vector of the flow over that edge -flow: Flow per flow edge in the order prescribed by flow_dict -flow_prev: The flow vector of the previous timestep, used for integration -flow_integrated: Flow integrated over time, used for mean flow computation - over saveat intervals saveat: The time interval between saves of output data (storage, flow, ...) """ const ModelGraph = MetaGraph{ @@ -778,17 +805,13 @@ const ModelGraph = MetaGraph{ node_ids::Dict{Int32, Set{NodeID}}, edges_source::Dict{Int32, Set{EdgeMetadata}}, flow_edges::Vector{EdgeMetadata}, - flow_dict::Dict{Tuple{NodeID, NodeID}, Int}, - flow::Cache, - flow_prev::Vector{Float64}, - flow_integrated::Vector{Float64}, saveat::Float64, }, MetaGraphsNext.var"#11#13", Float64, } -@kwdef struct Parameters{C1, C2, V1, V2} +@kwdef struct Parameters{C1, C2, C3, C4, V1, V2} starttime::DateTime graph::ModelGraph allocation::Allocation @@ -808,5 +831,17 @@ const ModelGraph = MetaGraph{ level_demand::LevelDemand flow_demand::FlowDemand subgrid::Subgrid + # Per state the in- and outflow edges associated with that state (if theu exist) + state_inflow_edge::C3 = ComponentVector() + state_outflow_edge::C4 = ComponentVector() all_nodes_active::Base.RefValue{Bool} = Ref(false) + tprev::Base.RefValue{Float64} = Ref(0.0) + # Water balance tolerances + water_balance_abstol::Float64 + water_balance_reltol::Float64 +end + +# To opt-out of type checking for ForwardDiff +function DiffEqBase.anyeltypedual(::Parameters, ::Type{Val{counter}}) where {counter} + Any end diff --git a/core/src/read.jl b/core/src/read.jl index e64e075b9..f57689956 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -571,6 +571,7 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin # both static and time are optional, but we need fallback defaults static = load_structvector(db, config, BasinStaticV1) time = load_structvector(db, config, BasinTimeV1) + state = load_structvector(db, config, BasinStateV1) set_static_value!(table, node_id, static) set_current_value!(table, node_id, time, config.starttime) @@ -578,15 +579,12 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin vertical_flux_from_input = ComponentVector(; precipitation, potential_evaporation, drainage, infiltration) - vertical_flux = cache(4 * n) - vertical_flux_prev = ComponentVector(; - precipitation = copy(precipitation), - evaporation, - drainage = copy(drainage), - infiltration = copy(infiltration), + vertical_flux_bmi = ComponentVector(; + precipitation = zero(precipitation), + evaporation = zero(evaporation), + drainage = zero(drainage), + infiltration = zero(infiltration), ) - vertical_flux_integrated = zero(vertical_flux_prev) - vertical_flux_bmi = zero(vertical_flux_prev) demand = zeros(length(node_id)) @@ -640,14 +638,11 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin error("Errors encountered when parsing Basin concentration data.") end - return Basin(; + basin = Basin(; node_id, inflow_ids = [collect(inflow_ids(graph, id)) for id in node_id], outflow_ids = [collect(outflow_ids(graph, id)) for id in node_id], vertical_flux_from_input, - vertical_flux, - vertical_flux_prev, - vertical_flux_integrated, vertical_flux_bmi, current_level, current_area, @@ -657,6 +652,12 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin time, concentration_external, ) + + storage0 = get_storages_from_levels(basin, state.level) + @assert length(storage0) == n "Basin / state length differs from number of Basins" + basin.storage0 .= storage0 + basin.storage_prev_saveat .= storage0 + return basin end """ @@ -697,7 +698,7 @@ function CompoundVariable( end function parse_variables_and_conditions(compound_variable, condition, ids, db, graph) - placeholder_vector = graph[].flow + placeholder_vector = cache(1) compound_variables = Vector{CompoundVariable}[] errors = false @@ -823,15 +824,8 @@ function continuous_control_functions(db, config, ids) return functions, controlled_variables, errors end -function continuous_control_compound_variables( - db::DB, - config::Config, - ids, - graph::MetaGraph, -) - # This is a vector that is known to have a DiffCache if automatic differentiation - # is used. Therefore this vector is used as a placeholder with the correct type - placeholder_vector = graph[].flow +function continuous_control_compound_variables(db::DB, config::Config, ids) + placeholder_vector = cache(1) data = load_structvector(db, config, ContinuousControlVariableV1) compound_variables = CompoundVariable[] @@ -860,7 +854,7 @@ function ContinuousControl(db::DB, config::Config, graph::MetaGraph)::Continuous # Avoid using `function` as a variable name as that is recognized as a keyword func, controlled_variable, errors = continuous_control_functions(db, config, ids) - compound_variable = continuous_control_compound_variables(db, config, ids, graph) + compound_variable = continuous_control_compound_variables(db, config, ids) # References to the controlled parameters, filled in later when they are known target_refs = PreallocationRef[] @@ -1248,7 +1242,6 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation end function Parameters(db::DB, config::Config)::Parameters - n_states = length(get_ids(db, "Basin")) + length(get_ids(db, "PidControl")) graph = create_graph(db, config) allocation = Allocation(db, config, graph) @@ -1298,6 +1291,8 @@ function Parameters(db::DB, config::Config)::Parameters level_demand, flow_demand, subgrid, + config.solver.water_balance_abstol, + config.solver.water_balance_reltol, ) collect_control_mappings!(p) diff --git a/core/src/solve.jl b/core/src/solve.jl index 7b1f2f02f..398adf5c6 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -7,13 +7,15 @@ function water_balance!( p::Parameters, t::Number, )::Nothing - (; graph, basin, pid_control) = p + (; basin, pid_control) = p du .= 0.0 - graph[].flow[parent(du)] .= 0.0 # Ensures current_* vectors are current - set_current_basin_properties!(basin, u, du) + set_current_basin_properties!(du, u, p, t) + + current_storage = basin.current_storage[parent(du)] + current_level = basin.current_level[parent(du)] # Notes on the ordering of these formulations: # - Continuous control can depend on flows (which are not continuously controlled themselves), @@ -22,10 +24,10 @@ function water_balance!( # because of the error derivative term. # Basin forcings - formulate_basins!(du, basin) + update_vertical_flux!(basin, du) # Formulate intermediate flows (non continuously controlled) - formulate_flows!(du, u, p, t) + formulate_flows!(du, p, t, current_storage, current_level) # Compute continuous control formulate_continuous_control!(du, p, t) @@ -33,38 +35,29 @@ function water_balance!( # Formulate intermediate flows (controlled by ContinuousControl) formulate_flows!( du, - u, p, - t; + t, + current_storage, + current_level; continuous_control_type = ContinuousControlType.Continuous, ) - # Formulate du (all) - formulate_du!(du, graph, u) - # Compute PID control formulate_pid_control!(u, du, pid_control, p, t) # Formulate intermediate flow (controlled by PID control) - formulate_flows!(du, u, p, t; continuous_control_type = ContinuousControlType.PID) - - # Formulate du (controlled by PidControl) - formulate_du_pid_controlled!(du, graph, pid_control) - - # https://github.com/Deltares/Ribasim/issues/1705#issuecomment-2283293974 - stop_declining_negative_storage!(du, u) + formulate_flows!( + du, + p, + t, + current_storage, + current_level; + continuous_control_type = ContinuousControlType.PID, + ) return nothing end -function stop_declining_negative_storage!(du, u) - for (i, s) in enumerate(u.storage) - if s < 0 - du.storage[i] = max(du.storage[i], 0.0) - end - end -end - function formulate_continuous_control!(du, p, t)::Nothing (; compound_variable, target_ref, func) = p.continuous_control @@ -75,63 +68,192 @@ function formulate_continuous_control!(du, p, t)::Nothing return nothing end +""" +Compute the storages, levels and areas of all Basins given the +state u and the time t. +""" function set_current_basin_properties!( - basin::Basin, - u::AbstractVector, - du::AbstractVector, + du::ComponentVector, + u::ComponentVector, + p::Parameters, + t::Number, )::Nothing - (; current_level, current_area) = basin + (; basin) = p + (; + current_storage, + current_level, + current_area, + current_cumulative_precipitation, + current_cumulative_drainage, + cumulative_precipitation, + cumulative_drainage, + vertical_flux_from_input, + ) = basin + current_storage = current_storage[parent(du)] current_level = current_level[parent(du)] current_area = current_area[parent(du)] + current_cumulative_precipitation = current_cumulative_precipitation[parent(du)] + current_cumulative_drainage = current_cumulative_drainage[parent(du)] + + # The exactly integrated precipitation and drainage up to the t of this water_balance call + dt = t - p.tprev[] + for node_id in basin.node_id + fixed_area = basin_areas(basin, node_id.idx)[end] + current_cumulative_precipitation[node_id.idx] = + cumulative_precipitation[node_id.idx] + + fixed_area * vertical_flux_from_input.precipitation[node_id.idx] * dt + end + @. current_cumulative_drainage = + cumulative_drainage + dt * vertical_flux_from_input.drainage + + formulate_storages!(current_storage, du, u, p, t) - for (i, s) in enumerate(u.storage) + for (i, s) in enumerate(current_storage) current_level[i] = get_level_from_storage(basin, i, s) current_area[i] = basin.level_to_area[i](current_level[i]) end end +function formulate_storages!( + current_storage::AbstractVector, + du::ComponentVector, + u::ComponentVector, + p::Parameters, + t::Number, +)::Nothing + (; + basin, + flow_boundary, + tabulated_rating_curve, + pump, + outlet, + linear_resistance, + manning_resistance, + user_demand, + tprev, + ) = p + # Current storage: initial conditdion + + # total inflows and outflows since the start + # of the simulation + current_storage .= basin.storage0 + formulate_storage!(current_storage, basin, du, u) + formulate_storage!(current_storage, tprev[], t, flow_boundary) + formulate_storage!(current_storage, t, u.tabulated_rating_curve, tabulated_rating_curve) + formulate_storage!(current_storage, t, u.pump, pump) + formulate_storage!(current_storage, t, u.outlet, outlet) + formulate_storage!(current_storage, t, u.linear_resistance, linear_resistance) + formulate_storage!(current_storage, t, u.manning_resistance, manning_resistance) + formulate_storage!( + current_storage, + t, + u.user_demand_inflow, + user_demand; + edge_volume_out = u.user_demand_outflow, + ) + return nothing +end + +""" +The storage contributions of the forcings that are part of the state. +""" +function formulate_storage!( + current_storage::AbstractVector, + basin::Basin, + du::ComponentVector, + u::ComponentVector, +) + (; current_cumulative_precipitation, current_cumulative_drainage) = basin + + current_storage .-= u.evaporation + current_storage .-= u.infiltration + + current_cumulative_precipitation = current_cumulative_precipitation[parent(du)] + current_cumulative_drainage = current_cumulative_drainage[parent(du)] + current_storage .+= current_cumulative_precipitation + current_storage .+= current_cumulative_drainage +end + +""" +Formulate storage contributions of nodes. +""" +function formulate_storage!( + current_storage::AbstractVector, + t::Number, + edge_volume_in::AbstractVector, + node::AbstractParameterNode; + edge_volume_out = nothing, +) + edge_volume_out = isnothing(edge_volume_out) ? edge_volume_in : edge_volume_out + + for (volume_in, volume_out, inflow_edge, outflow_edge) in + zip(edge_volume_in, edge_volume_out, node.inflow_edge, node.outflow_edge) + inflow_id = inflow_edge.edge[1] + if inflow_id.type == NodeType.Basin + current_storage[inflow_id.idx] -= volume_in + end + + outflow_id = outflow_edge.edge[2] + if outflow_id.type == NodeType.Basin + current_storage[outflow_id.idx] += volume_out + end + end +end + +""" +Formulate storage contributions of flow boundaries. +""" +function formulate_storage!( + current_storage::AbstractVector, + tprev::Number, + t::Number, + flow_boundary::FlowBoundary, +) + for (flow_rate, outflow_edges, active, cumulative_flow) in zip( + flow_boundary.flow_rate, + flow_boundary.outflow_edges, + flow_boundary.active, + flow_boundary.cumulative_flow, + ) + volume = cumulative_flow + if active + volume += integral(flow_rate, tprev, t) + end + for outflow_edge in outflow_edges + outflow_id = outflow_edge.edge[2] + if outflow_id.type == NodeType.Basin + current_storage[outflow_id.idx] += volume + end + end + end +end + """ Smoothly let the evaporation flux go to 0 when at small water depths Currently at less than 0.1 m. """ function update_vertical_flux!(basin::Basin, du::AbstractVector)::Nothing - (; current_level, current_area, vertical_flux_from_input, vertical_flux) = basin + (; vertical_flux_from_input, current_level, current_area) = basin current_level = current_level[parent(du)] current_area = current_area[parent(du)] - vertical_flux = wrap_forcing(vertical_flux[parent(du)]) for id in basin.node_id level = current_level[id.idx] area = current_area[id.idx] bottom = basin_levels(basin, id.idx)[1] - fixed_area = basin_areas(basin, id.idx)[end] depth = max(level - bottom, 0.0) factor = reduction_factor(depth, 0.1) - precipitation = fixed_area * vertical_flux_from_input.precipitation[id.idx] evaporation = area * factor * vertical_flux_from_input.potential_evaporation[id.idx] - drainage = vertical_flux_from_input.drainage[id.idx] infiltration = factor * vertical_flux_from_input.infiltration[id.idx] - vertical_flux.precipitation[id.idx] = precipitation - vertical_flux.evaporation[id.idx] = evaporation - vertical_flux.drainage[id.idx] = drainage - vertical_flux.infiltration[id.idx] = infiltration + du.evaporation[id.idx] = evaporation + du.infiltration[id.idx] = infiltration end return nothing end -function formulate_basins!(du::AbstractVector, basin::Basin)::Nothing - update_vertical_flux!(basin, du) - for id in basin.node_id - # add all vertical fluxes that enter the Basin - du.storage[id.idx] += get_influx(basin, id.idx) - end - return nothing -end - function set_error!(pid_control::PidControl, p::Parameters, du::ComponentVector, t::Number) (; basin) = p (; listen_node_id, target, error) = pid_control @@ -173,7 +295,7 @@ function formulate_pid_control!( listened_node_id = listen_node_id[i] - flow_rate = 0.0 + flow_rate = zero(eltype(du)) K_p = pid_control.proportional[i](t) K_i = pid_control.integral[i](t) @@ -198,11 +320,11 @@ function formulate_pid_control!( if !iszero(K_d) dlevel_demand = derivative(target[i], t) - du_listened_basin_old = du.storage[listened_node_id.idx] + dstorage_listened_basin_old = formulate_dstorage(du, p, t, listened_node_id) # The expression below is the solution to an implicit equation for - # du_listened_basin. This equation results from the fact that if the derivative + # dstorage_listened_basin. This equation results from the fact that if the derivative # term in the PID controller is used, the controlled pump flow rate depends on itself. - flow_rate += K_d * (dlevel_demand - du_listened_basin_old / area) / D + flow_rate += K_d * (dlevel_demand - dstorage_listened_basin_old / area) / D end # Set flow_rate @@ -211,17 +333,42 @@ function formulate_pid_control!( return nothing end +""" +Formulate the time derivative of the storage in a single Basin. +""" +function formulate_dstorage(du::ComponentVector, p::Parameters, t::Number, node_id::NodeID) + (; basin) = p + (; inflow_ids, outflow_ids, vertical_flux_from_input) = basin + @assert node_id.type == NodeType.Basin + dstorage = 0.0 + for inflow_id in inflow_ids[node_id.idx] + dstorage += get_flow(du, p, t, (inflow_id, node_id)) + end + for outflow_id in outflow_ids[node_id.idx] + dstorage -= get_flow(du, p, t, (node_id, outflow_id)) + end + + fixed_area = basin_areas(basin, node_id.idx)[end] + dstorage += fixed_area * vertical_flux_from_input.precipitation[node_id.idx] + dstorage += vertical_flux_from_input.drainage[node_id.idx] + dstorage -= du.evaporation[node_id.idx] + dstorage -= du.infiltration[node_id.idx] + + dstorage +end + function formulate_flow!( + du::ComponentVector, user_demand::UserDemand, - du::AbstractVector, - u::AbstractVector, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, )::Nothing - (; graph, allocation) = p + (; allocation) = p all_nodes_active = p.all_nodes_active[] for ( - node_id, + id, inflow_edge, outflow_edge, active, @@ -267,20 +414,17 @@ function formulate_flow!( # Smoothly let abstraction go to 0 as the source basin dries out inflow_id = inflow_edge.edge[1] - factor_basin = low_storage_factor(u.storage, inflow_id, 10.0) + factor_basin = low_storage_factor(current_storage, inflow_id, 10.0) q *= factor_basin # Smoothly let abstraction go to 0 as the source basin # level reaches its minimum level - source_level = get_level(p, inflow_id, t, du) + source_level = get_level(p, inflow_id, t, current_level) Δsource_level = source_level - min_level factor_level = reduction_factor(Δsource_level, 0.1) q *= factor_level - - set_flow!(graph, inflow_edge, q, du) - - # Return flow is immediate - set_flow!(graph, outflow_edge, q * return_factor(t), du) + du.user_demand_inflow[id.idx] = q + du.user_demand_outflow[id.idx] = q * return_factor(t) end return nothing end @@ -289,13 +433,13 @@ end Directed graph: outflow is positive! """ function formulate_flow!( + du::ComponentVector, linear_resistance::LinearResistance, - du::AbstractVector, - u::AbstractVector, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, )::Nothing - (; graph) = p all_nodes_active = p.all_nodes_active[] (; node_id, active, resistance, max_flow_rate) = linear_resistance for id in node_id @@ -306,14 +450,18 @@ function formulate_flow!( outflow_id = outflow_edge.edge[2] if (active[id.idx] || all_nodes_active) - h_a = get_level(p, inflow_id, t, du) - h_b = get_level(p, outflow_id, t, du) + h_a = get_level(p, inflow_id, t, current_level) + h_b = get_level(p, outflow_id, t, current_level) q_unlimited = (h_a - h_b) / resistance[id.idx] q = clamp(q_unlimited, -max_flow_rate[id.idx], max_flow_rate[id.idx]) - q *= low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, 10.0) - - set_flow!(graph, inflow_edge, q, du) - set_flow!(graph, outflow_edge, q, du) + q *= low_storage_factor_resistance_node( + current_storage, + q, + inflow_id, + outflow_id, + 10.0, + ) + du.linear_resistance[id.idx] = q end end return nothing @@ -323,13 +471,14 @@ end Directed graph: outflow is positive! """ function formulate_flow!( - tabulated_rating_curve::TabulatedRatingCurve, du::AbstractVector, - u::AbstractVector, + tabulated_rating_curve::TabulatedRatingCurve, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, )::Nothing - (; graph) = p + (; basin) = p all_nodes_active = p.all_nodes_active[] (; node_id, active, table) = tabulated_rating_curve @@ -340,12 +489,12 @@ function formulate_flow!( outflow_id = outflow_edge.edge[2] max_downstream_level = tabulated_rating_curve.max_downstream_level[id.idx] - h_a = get_level(p, inflow_id, t, du) - h_b = get_level(p, outflow_id, t, du) + h_a = get_level(p, inflow_id, t, current_level) + h_b = get_level(p, outflow_id, t, current_level) Δh = h_a - h_b if active[id.idx] || all_nodes_active - factor = low_storage_factor(u.storage, inflow_id, 10.0) + factor = low_storage_factor(current_storage, inflow_id, 10.0) q = factor * table[id.idx](h_a) q *= reduction_factor(Δh, 0.02) q *= reduction_factor(max_downstream_level - h_b, 0.02) @@ -353,8 +502,7 @@ function formulate_flow!( q = 0.0 end - set_flow!(graph, inflow_edge, q, du) - set_flow!(graph, outflow_edge, q, du) + du.tabulated_rating_curve[id.idx] = q end return nothing end @@ -399,13 +547,13 @@ hydraulic radius. This ensures that a basin can receive water after it has gone dry. """ function formulate_flow!( - manning_resistance::ManningResistance, du::AbstractVector, - u::AbstractVector, + manning_resistance::ManningResistance, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, )::Nothing - (; graph) = p (; node_id, active, @@ -428,8 +576,8 @@ function formulate_flow!( continue end - h_a = get_level(p, inflow_id, t, du) - h_b = get_level(p, outflow_id, t, du) + h_a = get_level(p, inflow_id, t, current_level) + h_b = get_level(p, outflow_id, t, current_level) bottom_a = upstream_bottom[id.idx] bottom_b = downstream_bottom[id.idx] @@ -457,50 +605,30 @@ function formulate_flow!( Δh = h_a - h_b q = A / n * ∛(R_h^2) * relaxed_root(Δh / L, 1e-3) - q *= low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, 10.0) - - set_flow!(graph, inflow_edge, q, du) - set_flow!(graph, outflow_edge, q, du) + q *= low_storage_factor_resistance_node( + current_storage, + q, + inflow_id, + outflow_id, + 10.0, + ) + du.manning_resistance[id.idx] = q end return nothing end function formulate_flow!( - flow_boundary::FlowBoundary, du::AbstractVector, - u::AbstractVector, - p::Parameters, - t::Number, -)::Nothing - (; graph, all_nodes_active) = p - all_nodes_active = p.all_nodes_active[] - (; node_id, active, flow_rate, outflow_edges) = flow_boundary - - for id in node_id - if active[id.idx] || all_nodes_active - rate = flow_rate[id.idx](t) - for outflow_edge in outflow_edges[id.idx] - - # Adding water is always possible - set_flow!(graph, outflow_edge, rate, du) - end - end - end -end - -function formulate_flow!( pump::Pump, - du::AbstractVector, - u::AbstractVector, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, continuous_control_type_::ContinuousControlType.T, )::Nothing - (; graph) = p all_nodes_active = p.all_nodes_active[] - for ( - node_id, + id, inflow_edge, outflow_edge, active, @@ -529,36 +657,33 @@ function formulate_flow!( inflow_id = inflow_edge.edge[1] outflow_id = outflow_edge.edge[2] - src_level = get_level(p, inflow_id, t, du) - dst_level = get_level(p, outflow_id, t, du) + src_level = get_level(p, inflow_id, t, current_level) + dst_level = get_level(p, outflow_id, t, current_level) - factor = low_storage_factor(u.storage, inflow_id, 10.0) + factor = low_storage_factor(current_storage, inflow_id, 10.0) q = flow_rate * factor q *= reduction_factor(src_level - min_upstream_level, 0.02) q *= reduction_factor(max_downstream_level - dst_level, 0.02) q = clamp(q, min_flow_rate, max_flow_rate) - - set_flow!(graph, inflow_edge, q, du) - set_flow!(graph, outflow_edge, q, du) + du.pump[id.idx] = q end return nothing end function formulate_flow!( - outlet::Outlet, du::AbstractVector, - u::AbstractVector, + outlet::Outlet, p::Parameters, t::Number, + current_storage::Vector, + current_level::Vector, continuous_control_type_::ContinuousControlType.T, )::Nothing - (; graph) = p all_nodes_active = p.all_nodes_active[] - for ( - node_id, + id, inflow_edge, outflow_edge, active, @@ -587,11 +712,11 @@ function formulate_flow!( inflow_id = inflow_edge.edge[1] outflow_id = outflow_edge.edge[2] - src_level = get_level(p, inflow_id, t, du) - dst_level = get_level(p, outflow_id, t, du) + src_level = get_level(p, inflow_id, t, current_level) + dst_level = get_level(p, outflow_id, t, current_level) q = flow_rate - q *= low_storage_factor(u.storage, inflow_id, 10.0) + q *= low_storage_factor(current_storage, inflow_id, 10.0) # No flow of outlet if source level is lower than target level Δlevel = src_level - dst_level @@ -600,73 +725,43 @@ function formulate_flow!( q *= reduction_factor(max_downstream_level - dst_level, 0.02) q = clamp(q, min_flow_rate, max_flow_rate) - - set_flow!(graph, inflow_edge, q, du) - set_flow!(graph, outflow_edge, q, du) - end - return nothing -end - -function formulate_du!(du::ComponentVector, graph::MetaGraph, u::AbstractVector)::Nothing - # loop over basins - # subtract all outgoing flows - # add all ingoing flows - for edge_metadata in values(graph[].flow_edges) - from_id, to_id = edge_metadata.edge - - if from_id.type == NodeType.Basin - q = get_flow(graph, edge_metadata, du) - du[from_id.idx] -= q - elseif to_id.type == NodeType.Basin - q = get_flow(graph, edge_metadata, du) - du[to_id.idx] += q - end - end - return nothing -end - -function formulate_du_pid_controlled!( - du::ComponentVector, - graph::MetaGraph, - pid_control::PidControl, -)::Nothing - for id in pid_control.controlled_basins - du[id.idx] = zero(eltype(du)) - for id_in in inflow_ids(graph, id) - du[id.idx] += get_flow(graph, id_in, id, du) - end - for id_out in outflow_ids(graph, id) - du[id.idx] -= get_flow(graph, id, id_out, du) - end + du.outlet[id.idx] = q end return nothing end function formulate_flows!( du::AbstractVector, - u::AbstractVector, p::Parameters, - t::Number; + t::Number, + current_storage::Vector, + current_level::Vector; continuous_control_type::ContinuousControlType.T = ContinuousControlType.None, )::Nothing (; linear_resistance, manning_resistance, tabulated_rating_curve, - flow_boundary, pump, outlet, user_demand, ) = p - formulate_flow!(pump, du, u, p, t, continuous_control_type) - formulate_flow!(outlet, du, u, p, t, continuous_control_type) + formulate_flow!(du, pump, p, t, current_storage, current_level, continuous_control_type) + formulate_flow!( + du, + outlet, + p, + t, + current_storage, + current_level, + continuous_control_type, + ) if continuous_control_type == ContinuousControlType.None - formulate_flow!(linear_resistance, du, u, p, t) - formulate_flow!(manning_resistance, du, u, p, t) - formulate_flow!(tabulated_rating_curve, du, u, p, t) - formulate_flow!(flow_boundary, du, u, p, t) - formulate_flow!(user_demand, du, u, p, t) + formulate_flow!(du, linear_resistance, p, t, current_storage, current_level) + formulate_flow!(du, manning_resistance, p, t, current_storage, current_level) + formulate_flow!(du, tabulated_rating_curve, p, t, current_storage, current_level) + formulate_flow!(du, user_demand, p, t, current_storage, current_level) end end diff --git a/core/src/util.jl b/core/src/util.jl index 610ce1b09..4bde59c6d 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -239,9 +239,8 @@ Get the current water level of a node ID. The ID can belong to either a Basin or a LevelBoundary. du: tells ForwardDiff whether this call is for differentiation or not """ -function get_level(p::Parameters, node_id::NodeID, t::Number, du)::Number +function get_level(p::Parameters, node_id::NodeID, t::Number, current_level::Vector)::Number if node_id.type == NodeType.Basin - current_level = p.basin.current_level[parent(du)] current_level[node_id.idx] elseif node_id.type == NodeType.LevelBoundary p.level_boundary.level[node_id.idx](t) @@ -356,7 +355,10 @@ function Base.getindex(fv::FlatVector, i::Int) end "Construct a FlatVector from one of the fields of SavedFlow." -FlatVector(saveval::Vector{SavedFlow}, sym::Symbol) = FlatVector(getfield.(saveval, sym)) +function FlatVector(saveval::Vector{<:SavedFlow}, sym::Symbol) + v = isempty(saveval) ? Vector{Float64}[] : getfield.(saveval, sym) + FlatVector(v) +end """ Function that goes smoothly from 0 to 1 in the interval [0,threshold], @@ -390,11 +392,17 @@ end For resistance nodes, give a reduction factor based on the upstream node as defined by the flow direction. """ -function low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, threshold) +function low_storage_factor_resistance_node( + current_storage, + q, + inflow_id, + outflow_id, + threshold, +) if q > 0 - low_storage_factor(u.storage, inflow_id, threshold) + low_storage_factor(current_storage, inflow_id, threshold) else - low_storage_factor(u.storage, outflow_id, threshold) + low_storage_factor(current_storage, outflow_id, threshold) end end @@ -565,31 +573,6 @@ function get_Δt(integrator)::Float64 end end -function get_influx(basin::Basin, node_id::NodeID)::Float64 - if node_id.type !== NodeType.Basin - error("Sum of vertical fluxes requested for non-basin $node_id.") - end - return get_influx(basin, node_id.idx) -end - -function get_influx(basin::Basin, basin_idx::Int; prev::Bool = false)::Float64 - (; node_id, vertical_flux, vertical_flux_prev, vertical_flux_from_input) = basin - influx = if prev - vertical_flux_prev.precipitation[basin_idx] - - vertical_flux_prev.evaporation[basin_idx] + - vertical_flux_prev.drainage[basin_idx] - - vertical_flux_prev.infiltration[basin_idx] - else - n = length(node_id) - vertical_flux = vertical_flux[parent(vertical_flux_from_input)] - vertical_flux[basin_idx] - # precipitation - vertical_flux[n + basin_idx] + # evaporation - vertical_flux[2n + basin_idx] - # drainage - vertical_flux[3n + basin_idx] # infiltration - end - return influx -end - inflow_edge(graph, node_id)::EdgeMetadata = graph[inflow_id(graph, node_id), node_id] outflow_edge(graph, node_id)::EdgeMetadata = graph[node_id, outflow_id(graph, node_id)] outflow_edges(graph, node_id)::Vector{EdgeMetadata} = @@ -613,9 +596,9 @@ function set_initial_allocation_mean_flows!(integrator)::Nothing for edge in keys(mean_input_flows) if edge[1] == edge[2] - q = get_influx(p.basin, edge[1]) + q = get_influx(du, edge[1], p) else - q = get_flow(graph, edge..., du) + q = get_flow(du, p, t, edge) end # Multiply by Δt_allocation as averaging divides by this factor # in update_allocation! @@ -628,7 +611,7 @@ function set_initial_allocation_mean_flows!(integrator)::Nothing if edge[1] == edge[2] mean_realized_flows[edge] = -u[edge[1].idx] else - q = get_flow(graph, edge..., du) + q = get_flow(du, p, t, edge) mean_realized_flows[edge] = q * Δt_allocation end end @@ -662,6 +645,9 @@ function get_variable_ref( (; basin, graph) = p errors = false + # Only built here because it is needed to obtain indices + u = build_state_vector(p) + ref = if node_id.type == NodeType.Basin && variable == "level" PreallocationRef(basin.current_level, node_id.idx) elseif variable == "flow_rate" && node_id.type != NodeType.FlowBoundary @@ -671,9 +657,9 @@ function get_variable_ref( @error "Cannot listen to flow_rate of $node_id, the node type must be one of $conservative_node_types" Ref(Float64[], 0) else - id_in = inflow_id(graph, node_id) - (; flow_idx) = graph[id_in, node_id] - PreallocationRef(graph[].flow, flow_idx) + # Index in the state vector (inflow) + flow_idx = get_state_index(node_id, u) + PreallocationRef(cache(1), flow_idx; from_du = true) end else node = getfield(p, snake_case(Symbol(node_id.type))) @@ -681,7 +667,7 @@ function get_variable_ref( end else # Placeholder to obtain correct type - PreallocationRef(graph[].flow, 0) + PreallocationRef(cache(1), 0) end return ref, errors end @@ -848,37 +834,6 @@ function basin_areas(basin::Basin, state_idx::Int) return basin.level_to_area[state_idx].u end -""" -Just before setting a timestep, call water_balance! again -to get a correct value of the flows for integrating -""" -function set_previous_flows!(integrator) - (; p, u, t) = integrator - (; flow, flow_prev) = p.graph[] - (; vertical_flux, vertical_flux_prev) = p.basin - du = get_du(integrator) - water_balance!(du, u, p, t) - - flow = flow[parent(u)] - vertical_flux = vertical_flux[parent(u)] - copyto!(flow_prev, flow) - copyto!(vertical_flux_prev, vertical_flux) -end - -""" -Split the single forcing vector into views of the components -precipitation, evaporation, drainage, infiltration -""" -function wrap_forcing(vector) - n = length(vector) ÷ 4 - (; - precipitation = view(vector, 1:n), - evaporation = view(vector, (n + 1):(2n)), - drainage = view(vector, (2n + 1):(3n)), - infiltration = view(vector, (3n + 1):(4n)), - ) -end - """ The function f(x) = sign(x)*√(|x|) where for |x| model.integrator.sol(t)[1] + itp_basin_2 = LinearInterpolation(storage, t) realized_numeric = diff(itp_basin_2.(df_basin_2.time)) / Δt_allocation @test all(isapprox.(realized_numeric, df_basin_2.realized[2:end], atol = 2e-4)) # Realized user demand flow_table = DataFrame(Ribasim.flow_table(model)) - flow_table_user_3 = flow_table[flow_table.edge_id .== 1, :] + flow_table_user_3 = flow_table[flow_table.edge_id .== 2, :] itp_user_3 = LinearInterpolation( flow_table_user_3.flow_rate, Ribasim.seconds_since.(flow_table_user_3.time, model.config.starttime), diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index d7cdd036e..6e0a51800 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -4,7 +4,7 @@ toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) @test BMI.get_time_units(model) == "s" - dt0 = 0.0001269439f0 + dt0 = 2.8280652f-5 @test BMI.get_time_step(model) ≈ dt0 atol = 5e-3 @test BMI.get_start_time(model) === 0.0 @test BMI.get_current_time(model) === 0.0 @@ -93,11 +93,12 @@ end slope = 1e-3 / year day = 86400.0 BMI.update_until(model, day) - @test realized ≈ [demand_start * day, demand_start * day + 0.5 * slope * day^2] + @test realized ≈ [demand_start * day, demand_start * day + 0.5 * slope * day^2] atol = + 1e-3 demand_later = 2e-3 demand[1] = demand_later BMI.update_until(model, 2day) - @test realized[1] == demand_start * day + demand_later * day + @test realized[1] ≈ demand_start * day + demand_later * day atol = 1e-3 end @testitem "vertical basin flux" begin diff --git a/core/test/control_test.jl b/core/test/control_test.jl index 937633171..3109cc1e7 100644 --- a/core/test/control_test.jl +++ b/core/test/control_test.jl @@ -1,5 +1,6 @@ @testitem "Pump discrete control" begin using Ribasim: NodeID + using OrdinaryDiffEqCore: get_du using Dates: DateTime import Arrow import Tables @@ -56,7 +57,7 @@ t_2_index = findfirst(>=(t_2), t) @test level[2, t_2_index] >= discrete_control.compound_variables[1][2].greater_than[1] - flow = graph[].flow[Float64[]] + flow = get_du(model.integrator)[(:linear_resistance, :pump)] @test all(iszero, flow) end @@ -247,7 +248,13 @@ end t_switch = Ribasim.datetime_since(record.time[2], p.starttime) flow_table = DataFrame(Ribasim.flow_table(model)) @test all(filter(:time => time -> time <= t_switch, flow_table).flow_rate .> 0) - @test all(filter(:time => time -> time > t_switch, flow_table).flow_rate .== 0) + @test all( + isapprox.( + filter(:time => time -> time > t_switch, flow_table).flow_rate, + 0; + atol = 1e-8, + ), + ) end @testitem "Outlet continuous control" begin diff --git a/core/test/docs.toml b/core/test/docs.toml index e16a9d775..1ab1721d3 100644 --- a/core/test/docs.toml +++ b/core/test/docs.toml @@ -33,6 +33,8 @@ dtmax = 0.0 # optional, default length of simulation force_dtmin = false # optional, default false abstol = 1e-6 # optional, default 1e-6 reltol = 1e-5 # optional, default 1e-5 +water_balance_abstol = 1e-6 # optional, default 1e-6 +water_balance_reltol = 1e-6 # optional, default 1e-6 maxiters = 1e9 # optional, default 1e9 sparse = true # optional, default true autodiff = true # optional, default true diff --git a/core/test/io_test.jl b/core/test/io_test.jl index 5bb70e92c..45d552a94 100644 --- a/core/test/io_test.jl +++ b/core/test/io_test.jl @@ -97,6 +97,7 @@ end # load a sorted table table = Ribasim.load_structvector(db, config, Ribasim.BasinTimeV1) + close(db) by = Ribasim.sort_by_function(table) @test by == Ribasim.sort_by_time_id # reverse it so it needs sorting @@ -153,9 +154,9 @@ end config = Ribasim.Config(toml_path) model = Ribasim.Model(config) - storage1_begin = copy(model.integrator.u.storage) + storage1_begin = copy(model.integrator.p.basin.current_storage[Float64[]]) solve!(model) - storage1_end = model.integrator.u.storage + storage1_end = model.integrator.p.basin.current_storage[Float64[]] @test storage1_begin != storage1_end # copy state results to input @@ -171,6 +172,6 @@ end end model = Ribasim.Model(toml_path) - storage2_begin = model.integrator.u.storage + storage2_begin = model.integrator.p.basin.current_storage[Float64[]] @test storage1_end ≈ storage2_begin end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index c0627f766..349546f50 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -123,22 +123,25 @@ end @testitem "bucket model" begin using SciMLBase: successful_retcode + using OrdinaryDiffEqCore: get_du toml_path = normpath(@__DIR__, "../../generated_testmodels/bucket/ribasim.toml") @test ispath(toml_path) model = Ribasim.run(toml_path) @test model isa Ribasim.Model - @test model.integrator.u.storage ≈ [1000] - vertical_flux = Ribasim.wrap_forcing(model.integrator.p.basin.vertical_flux[Float64[]]) - @test vertical_flux.precipitation == [0.0] - @test vertical_flux.evaporation == [0.0] - @test vertical_flux.drainage == [0.0] - @test vertical_flux.infiltration == [0.0] + (; basin) = model.integrator.p + @test basin.current_storage[Float64[]] ≈ [1000] + @test basin.vertical_flux_from_input.precipitation == [0.0] + @test basin.vertical_flux_from_input.drainage == [0.0] + du = get_du(model.integrator) + @test du.evaporation == [0.0] + @test du.infiltration == [0.0] @test successful_retcode(model) end @testitem "leaky bucket model" begin using SciMLBase: successful_retcode + using OrdinaryDiffEqCore: get_du import BasicModelInterface as BMI toml_path = normpath(@__DIR__, "../../generated_testmodels/leaky_bucket/ribasim.toml") @@ -146,12 +149,15 @@ end model = Ribasim.Model(toml_path) @test model isa Ribasim.Model - stor = model.integrator.u.storage - vertical_flux = Ribasim.wrap_forcing(model.integrator.p.basin.vertical_flux[Float64[]]) - prec = vertical_flux.precipitation - evap = vertical_flux.evaporation - drng = vertical_flux.drainage - infl = vertical_flux.infiltration + (; integrator) = model + du = get_du(integrator) + (; u, p, t) = integrator + Ribasim.water_balance!(du, u, p, t) + stor = integrator.p.basin.current_storage[parent(du)] + prec = p.basin.vertical_flux_from_input.precipitation + evap = du.evaporation + drng = p.basin.vertical_flux_from_input.drainage + infl = du.infiltration # The dynamic data has missings, but these are not set. @test prec == [0.0] @test evap == [0.0] @@ -198,8 +204,8 @@ end @test successful_retcode(model) @test allunique(Ribasim.tsaves(model)) - @test model.integrator.u ≈ Float32[803.7093, 803.68274, 495.241, 1318.3053] skip = - Sys.isapple() atol = 1.5 + @test model.integrator.p.basin.current_storage[Float64[]] ≈ + Float32[803.7093, 803.68274, 495.241, 1318.3053] skip = Sys.isapple() atol = 1.5 @test length(logger.logs) > 10 @test logger.logs[1].level == Debug @@ -223,6 +229,7 @@ end @testitem "basic transient model" begin using SciMLBase: successful_retcode + using OrdinaryDiffEqCore: get_du toml_path = normpath(@__DIR__, "../../generated_testmodels/basic_transient/ribasim.toml") @@ -231,13 +238,11 @@ end @test model isa Ribasim.Model @test successful_retcode(model) @test allunique(Ribasim.tsaves(model)) - precipitation = - Ribasim.wrap_forcing( - model.integrator.p.basin.vertical_flux[Float64[]], - ).precipitation + du = get_du(model.integrator) + precipitation = model.integrator.p.basin.vertical_flux_from_input.precipitation @test length(precipitation) == 4 - @test model.integrator.u ≈ Float32[698.22736, 698.2014, 421.20447, 1334.4354] atol = 2.0 skip = - Sys.isapple() + @test model.integrator.p.basin.current_storage[parent(du)] ≈ + Float32[697.30591, 697.2799, 419.19034, 1334.3859] atol = 2.0 skip = Sys.isapple() end @testitem "Allocation example model" begin @@ -278,7 +283,7 @@ end config = Ribasim.Config(toml_path; solver_algorithm = "Rodas5", solver_autodiff = true) time_ad = Ribasim.run(config) @test successful_retcode(time_ad) - @test time_ad.integrator.u ≈ sparse_ad.integrator.u atol = 4 + @test time_ad.integrator.u ≈ sparse_ad.integrator.u atol = 10 end @testitem "TabulatedRatingCurve model" begin @@ -290,7 +295,8 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test successful_retcode(model) - @test model.integrator.u ≈ Float32[368.31558, 365.68442] skip = Sys.isapple() + @test model.integrator.p.basin.current_storage[Float64[]] ≈ + Float32[368.31558, 365.68442] skip = Sys.isapple() # the highest level in the dynamic table is updated to 1.2 from the callback @test model.integrator.p.tabulated_rating_curve.table[end].t[end] == 1.2 end @@ -390,18 +396,25 @@ end using SciMLBase: successful_retcode using Dates using DataFrames: DataFrame + using Ribasim: formulate_storages! + import BasicModelInterface as BMI toml_path = normpath(@__DIR__, "../../generated_testmodels/user_demand/ribasim.toml") @test ispath(toml_path) - model = Ribasim.run(toml_path) - @test successful_retcode(model) + model = Ribasim.Model(toml_path) + + (; u, p, t, sol) = model.integrator + current_storage = p.basin.current_storage[Float64[]] - seconds_in_day = 86400.0 - @test only(model.integrator.sol(0seconds_in_day)) == 1000.0 + day = 86400.0 + + @test only(current_storage) ≈ 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 + BMI.update_until(model, 150day) + @test only(current_storage) ≈ 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 + BMI.update_until(model, 220day) + @test only(current_storage) ≈ 500 atol = 1 # Trasient return factor flow = DataFrame(Ribasim.flow_table(model)) @@ -421,6 +434,7 @@ end @testitem "ManningResistance" begin using PreallocationTools: get_tmp using SciMLBase: successful_retcode + using OrdinaryDiffEqCore: get_du using Ribasim: NodeID """ @@ -482,9 +496,9 @@ end model = Ribasim.run(toml_path) @test successful_retcode(model) - u = model.integrator.sol.u[end] - p = model.integrator.p - h_actual = p.basin.current_level[parent(u)][1:50] + du = get_du(model.integrator) + (; p, t) = model.integrator + h_actual = p.basin.current_level[parent(du)][1:50] x = collect(10.0:20.0:990.0) h_expected = standard_step_method(x, 5.0, 1.0, 0.04, h_actual[end], 1.0e-6) @@ -494,18 +508,13 @@ end # https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/latest/theoretical-basis-for-one-dimensional-and-two-dimensional-hydrodynamic-calculations/1d-steady-flow-water-surface-profiles/friction-loss-evaluation @test all(isapprox.(h_expected, h_actual; atol = 0.02)) # Test for conservation of mass, flow at the beginning == flow at the end - n_self_loops = length(p.graph[].flow_dict) - @test Ribasim.get_flow( - p.graph, - NodeID(:FlowBoundary, 1, p), - NodeID(:Basin, 2, p), - parent(u), - ) ≈ 5.0 atol = 0.001 skip = Sys.isapple() + @test Ribasim.get_flow(du, p, t, (NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))) ≈ + 5.0 atol = 0.001 skip = Sys.isapple() @test Ribasim.get_flow( - p.graph, - NodeID(:ManningResistance, 101, p), - NodeID(:Basin, 102, p), - parent(u), + du, + p, + t, + (NodeID(:ManningResistance, 101, p), NodeID(:Basin, 102, p)), ) ≈ 5.0 atol = 0.001 skip = Sys.isapple() end diff --git a/core/test/time_test.jl b/core/test/time_test.jl index 76ba54f12..14b535117 100644 --- a/core/test/time_test.jl +++ b/core/test/time_test.jl @@ -29,6 +29,7 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/basic_transient/ribasim.toml") @test ispath(toml_path) + config = Ribasim.Config(toml_path; solver_saveat = 0) model = Ribasim.run(toml_path) (; p) = model.integrator (; basin) = p @@ -50,9 +51,8 @@ end get_area(basin, idx, storage) 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 + # Mean areas over the timestep are computed as an approximation of + # how the area changes over the timestep, affecting evaporation time_table[!, "mean_area"] .= 0.0 n_basins = length(basin.node_id) n_times = length(unique(time_table.time)) - 1 @@ -69,7 +69,7 @@ end isapprox( basin_table.evaporation, time_table.mean_area .* time_table.potential_evaporation; - rtol = 1e-4, + rtol = 1e-2, ), ) @@ -96,11 +96,11 @@ end solver_algorithm = "Euler", ) model = Ribasim.Model(config) - starting_precipitation = Ribasim.wrap_forcing( - model.integrator.p.basin.vertical_flux[Float64[]], - ).precipitation[1] + (; basin) = model.integrator.p + starting_precipitation = + basin.vertical_flux_from_input.precipitation[1] * Ribasim.basin_areas(basin, 1)[end] BMI.update_until(model, saveat) - mean_precipitation = only(model.saved.vertical_flux.saveval).precipitation[1] + mean_precipitation = only(model.saved.flow.saveval).precipitation[1] # Given that precipitation stops after 15 of the 20 days @test mean_precipitation ≈ 3 / 4 * starting_precipitation end diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index a95a3e7e3..8d5b54198 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -161,7 +161,7 @@ end @testitem "Jacobian sparsity" begin import SQLite using ComponentArrays: ComponentVector - using SparseArrays: spzeros + using SparseArrays: sparse, findnz toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") @@ -170,17 +170,19 @@ end db = SQLite.DB(db_path) p = Ribasim.Parameters(db, cfg) + close(db) t0 = 0.0 - u0 = ComponentVector{Float64}(; - storage = zeros(length(p.basin.node_id)), - integral = Float64[], - ) + u0 = Ribasim.build_state_vector(p) du0 = copy(u0) jac_prototype = Ribasim.get_jac_prototype(du0, u0, p, t0) - jac_prototype_expected = spzeros(Bool, 4, 4) - jac_prototype_expected[1:2, 1:2] .= true - jac_prototype_expected[2:4, 2:4] .= true + # rows, cols, _ = findnz(jac_prototype) + #! format: off + rows_expected = [1, 2, 3, 6, 7, 9, 13, 1, 2, 3, 4, 6, 7, 9, 10, 13, 14, 1, 2, 3, 4, 5, 6, 7, 9, 11, 13, 15, 2, 3, 4, 5, 10, 11, 14, 15, 3, 4, 5, 11, 15, 1, 2, 3, 6, 7, 9, 13, 1, 2, 3, 6, 7, 8, 9, 12, 13, 7, 8, 12, 1, 2, 3, 6, 7, 9, 13, 2, 4, 10, 14, 3, 4, 5, 11, 15, 7, 8, 12, 1, 2, 3, 6, 7, 9, 13, 2, 4, 10, 14, 3, 4, 5, 11, 15] + cols_expected = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 15] + #! format: on + jac_prototype_expected = + sparse(rows_expected, cols_expected, true, size(jac_prototype)...) @test jac_prototype == jac_prototype_expected toml_path = normpath(@__DIR__, "../../generated_testmodels/pid_control/ribasim.toml") @@ -190,16 +192,17 @@ end db = SQLite.DB(db_path) p = Ribasim.Parameters(db, cfg) - u0 = ComponentVector{Float64}(; - storage = zeros(length(p.basin.node_id)), - integral = zeros(length(p.pid_control.node_id)), - ) + close(db) + u0 = Ribasim.build_state_vector(p) du0 = copy(u0) jac_prototype = Ribasim.get_jac_prototype(du0, u0, p, t0) - jac_prototype_expected = spzeros(Bool, 3, 3) - jac_prototype_expected[1, 1:3] .= true - jac_prototype_expected[2:3, 1] .= true + #! format: off + rows_expected = [1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2] + cols_expected = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 6] + #! format: on + jac_prototype_expected = + sparse(rows_expected, cols_expected, true, size(jac_prototype)...) @test jac_prototype == jac_prototype_expected end @@ -246,27 +249,12 @@ end node_id_1 = NodeID(:Basin, 5, 1) node_id_2 = NodeID(:Basin, 6, 2) - @test low_storage_factor_resistance_node( - (; storage = [3.0, 3.0]), - 1.0, - node_id_1, - node_id_2, - 2.0, - ) == 1.0 - @test low_storage_factor_resistance_node( - (; storage = [1.0, 3.0]), - 1.0, - node_id_1, - node_id_2, - 2.0, - ) == 0.5 - @test low_storage_factor_resistance_node( - (; storage = [1.0, 3.0]), - -1.0, - node_id_1, - node_id_2, - 2.0, - ) == 1.0 + @test low_storage_factor_resistance_node([3.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) == + 1.0 + @test low_storage_factor_resistance_node([1.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) == + 0.5 + @test low_storage_factor_resistance_node([1.0, 3.0], -1.0, node_id_1, node_id_2, 2.0) == + 1.0 end @testitem "constraints_from_nodes" begin diff --git a/core/test/validation_test.jl b/core/test/validation_test.jl index 1e8d73060..737d892f5 100644 --- a/core/test/validation_test.jl +++ b/core/test/validation_test.jl @@ -88,7 +88,6 @@ end function set_edge_metadata!(id_1, id_2, edge_type) graph[id_1, id_2] = EdgeMetadata(; id = 0, - flow_idx = 0, type = edge_type, subnetwork_id_source = 0, edge = (id_1, id_2), @@ -151,7 +150,6 @@ end function set_edge_metadata!(id_1, id_2, edge_type) graph[id_1, id_2] = EdgeMetadata(; id = 0, - flow_idx = 0, type = edge_type, subnetwork_id_source = 0, edge = (id_1, id_2), @@ -199,6 +197,7 @@ end db_path = Ribasim.input_path(cfg, cfg.database) db = SQLite.DB(db_path) p = Ribasim.Parameters(db, cfg) + close(db) logger = TestLogger() with_logger(logger) do @@ -281,6 +280,7 @@ end with_logger(logger) do @test !Ribasim.valid_edge_types(db) end + close(db) @test length(logger.logs) == 2 @test logger.logs[1].level == Error @@ -446,7 +446,9 @@ end "Warning: Convergence bottlenecks in descending order of severity:", output, ) - @test occursin("Basin #11 = ", output) + @test occursin("Pump #12 = ", output) + @test occursin("Pump #32 = ", output) + @test occursin("Pump #52 = ", output) end @testitem "Missing priority when allocation is active" begin diff --git a/docs/concept/allocation.qmd b/docs/concept/allocation.qmd index 482bc82e2..0901ba7d3 100644 --- a/docs/concept/allocation.qmd +++ b/docs/concept/allocation.qmd @@ -251,7 +251,6 @@ allocation_model = p.allocation.allocation_models[1] t = 0.0 priority_idx = 1 -Ribasim.set_flow!(p.graph, 1, 1.0, u) Ribasim.set_objective_priority!(allocation_model, p, u, t, priority_idx) Ribasim.set_initial_values!(allocation_model, p, u, t) diff --git a/docs/concept/equations.qmd b/docs/concept/equations.qmd index 36aff1b81..ed03e89c6 100644 --- a/docs/concept/equations.qmd +++ b/docs/concept/equations.qmd @@ -3,76 +3,95 @@ title: "Equations" --- # Formal model description -In this section we give a formal description of the problem that is solved by Ribasim. -The problem is of the form - +In this section we give a formal description of the problem that is solved by Ribasim. The problem is of the form $$ -\frac{\text{d}\mathbf{u}}{\text{d}t} = f(\mathbf{u},p(t),t),\quad t \in [t_0,t_\text{end}], + \frac{\text{d}\mathbf{u}}{\text{d}t} = f(\mathbf{u},p(t),t), \quad t \in [t_0, t_\text{end}], $$ -i.e. a system of coupled first order ordinary differential equations, with initial condition $\mathbf{u}(t_0)= \mathbf{u}_0$ and time dependent input data denoted by $p(t)$. - -The model is given by a directed graph, consisting of a set of nodes (or vertices) $V$ and edges $E$. -Let $V$ be the set of node IDs and let $E$ be the set of ordered tuples $(i,j)$ meaning that node $i$ is connected to node $j$. - -We can split the set of nodes into two subsets $V = B \cup N$, where $B$ is the set of basins and $N$ is the set of non-basins. -The basins have an associated storage state and the non-basins dictate how water flows to or from basins. +which is a system of coupled first order differential equations. -$\mathbf{u}(t)$ is given by all the states of the model, which are (currently) the storage of the basins and the integral terms of the PID controllers, the latter being explained in [PID equations](/reference/node/pid-control.qmd#equations). - -Given a single basin with node ID $i \in B$, the equation that dictates the change of its storage over time is given by +The model is given by a directed graph, consisting of a set of node IDs (vertices) $V$ and edges $E$, consisting of ordered pairs of node IDs. +We denote the subset of the nodes given by the Basins $B \subset V$, and the subset of nodes that prescribe flow $N \subset V$. +The states $\mathbf{u}$ of the model are given by cumulative flows since the start of the simulation as prescribed by the nodes $N$: $$ -\frac{\text{d}u_i}{\text{d}t} = -\sum_{(i',j') \in E | j' = i} Q_{i',j'} - \sum_{(i',j') \in E | i' = i} Q_{i',j'} + F_i(p,t). + u_n(t) = \int_{t_0}^t q_n\text{d}t' \quad \forall n \in N, +$$ +as well as by the Basin forcings: +$$ + u_b^\text{forcing}(t) = \int_{t_0}^t q_b^\text{forcing}\text{d}t' \quad \forall b \in B. $$ -Here $Q_{i,j}$ is the flow along an edge, where the graph direction dictates positive flow. -So the first term denotes flow towards the basin, the second one denotes flow away from the basin, and the third term denotes external forcing. -$F_i(p,t)$ is given by input data, and $Q_{i' ,j'}$ is determined by the type of nodes that connect to that edge. - -The various node and forcing types that the model can contain are explained in the section [Natural water balance terms](/concept/equations.qmd#natural-water-balance-terms). - -::: {.callout-note} -In general a model has more nodes than states, so in the Julia core there is a distinction between node indices and state indices. For simplicity these are treated as equal in the documentation when it comes to basins and their storage. -::: +Because of this definition, the initial conditions of all states are simple: +$$ + u_i(t_0) = 0 \quad \forall i. +$$ -## The Jacobian +From these cumulative flows, the storage in each Basin can be determined at each point in time: +$$ + S_b(t) = S_i(0) + S^\text{exact}(t) - u_b^\text{forcing}(t) + \sum_{n\;|\;(n,b)\in E} u(t) - \sum_{n\;|\;(b,n)\in E} u(t), +$$ -The Jacobian is a $n\times n$ matrix where $n$ is the number of states in the simulation. The Jacobian is computed either using finite difference methods or automatic differentiation. For more details on the computation of the Jacobian and how it is used in the solvers see [numerical considerations](numerics.qmd). +i.e. the storage is given by: +- the initial storage; +- plus the exactly integrated flows (more on that below); +- minus the cumulative outgoing forcings; +- plus the cumulative horizontal inflows; +- minus the cumulative horizontal outflows. -The entries of the Jacobian $J$ are given by +From these storages in combination with the Basin profiles the Basin levels $h$ are computed. +The relationship between the profile and the storage is given by $$ -J[i,j] = \frac{\partial f_j}{\partial u_i}, + S_b = \int_{h_0}^h A_b(\ell)\text{d}\ell, $$ -hence $J[i,j]$ quantifies how $f_j$, the derivative of state $j$ with respect to time, changes with a change in state $i$. If a node creates dependendies between basin storages (or other states), then this yields contributions to the Jacobian. If $j$ corresponds to a storage state, then +where $A_b$ is the linear interpolation of the area as a function of the level. +These levels are then inputs for determining the flows prescribed by the nodes $N$. From this relation it also follows that $$ -J[i,j] = \sum_{(i',j') \in E | j' = i} \frac{\partial Q_{i',j'}}{\partial u_i} - \sum_{(i',j') \in E | i' = i} \frac{\partial Q_{i',j'}}{\partial u_i}, + \frac{\text{d}h}{\text{d}t} = \frac{1}{A_b}, $$ +and so areas of zero are not allowed in the Basin profiles. -Most of these terms are always $0$, because a flow over an edge only depends on a small number of states. Therefore the matrix $J$ is very sparse. +## The PID control integral state -For many contributions to the Jacobian the derivative of the level $l(u)$ of a basin with respect to its storage $u$ is required. To get an expression for this, we first look at the storage as a function of the level: +There's one other type of state, which is not a cumulative flow but a cumulative error. +This is the error integral for PID control, further explained in [PID equations](/reference/node/pid-control.qmd#equations). -$$ -u(l) = \int_{l_0}^l A(\ell)d\ell. -$$ +## Exactly integrating flows to minimize the number of states + +The more states the problem has, the more time it takes to solve it. +Therefore we want to minimize the number of states. +Flows do not have to be states when they can be integrated over time exactly because they do not depend on the other states. +This is true for FlowBoundary nodes, and Basin precipitation (which uses a fixed basin area) and drainage. -From this we obtain $u'(l) = A(l)$ and thus +## The Jacobian + +The Jacobian is an $N \times N$ matrix where $N$ is the number of states in the simulation. +It is computed as part of implicit time stepping methods. +There are 2 different methods available for computing this matrix: finite difference or automatic differentiation. +For more details on the computation of the Jacobian and how it is used in the solvers see [numerical considerations](numerics.qmd). + +The entries of the Jacobian $J$ are given by $$ -\frac{\text{d}l}{\text{d}u} = \frac{1}{A(u)}. + J_{i,j} = \frac{\partial f_j}{\partial u_i}, $$ +i.e. $J_{i,j}$ quantifies how $f_j$. the time derivative of state $j$, changes with respect to changes in state $i$. Most of these entries are $0$, because flows in distant parts of the model do not depend on each other. + +## Why this formulation + +You might wonder why in the above explanation the states are given by the cumulative flows and not by the Basin storages, which is arguably conceptually simpler. +The reason is that we do not just want to model the storages in the Basins over time, but we also want accurate output of each individual flow, e.g. to model the spread of pollutants. -:::{.callout-note} -The presence of division by the basin area means that areas of size zero are not allowed. -::: +When the states are given by the storages, generally the individual flows can not accurately be computed from that as a post processing step, because there are more flows than storages. +Also, we can only compute flows at individual points in time explicitly, not over a whole interval. +When the states are given by the cumulative flows however, the output of the problem solve gives these flows directly, and from those the storage over time can be computed accurately. +Hence in short, the formulation above gives more information than a formulation with Basin storages as states. -# Numerical solution +## Numerical solution -Ribasim uses OrdinaryDiffEq.jl to provide a numerical solution to the water balance equations. +Ribasim uses [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl/) to provide a numerical solution to the water balance equations. Changes to forcings or parameters such as precipitation, but also the allocated water abstraction is managed through the use of callback functions [@callbacks]. In a coupled run, the exchanges with MODFLOW 6 are also managed via the use of a callback function. For more a more in-depth discussion of numerical computations see [Numerical considerations](numerics.qmd). diff --git a/docs/concept/modelconcept.qmd b/docs/concept/modelconcept.qmd index 96361f5f7..c9e46bea9 100644 --- a/docs/concept/modelconcept.qmd +++ b/docs/concept/modelconcept.qmd @@ -26,6 +26,10 @@ $$ \frac{\mathrm{d}S}{\mathrm{d}t} = P + ET + Q_{rest} $$ +We don't use these equations directly. +Rather, we use an equivalent formulation where solve for the cumulative flows instead of the Basin storages. +For more details on this see [Equations](equations.qmd). + ## Time The water balance equation can be applied on many timescales; years, weeks, days or hours. diff --git a/docs/concept/numerics.qmd b/docs/concept/numerics.qmd index 9a6fcee5a..35e887204 100644 --- a/docs/concept/numerics.qmd +++ b/docs/concept/numerics.qmd @@ -10,13 +10,14 @@ $$ \end{cases}, $$ {#eq-prob} -where $\mathbf{f}$ denotes `water_balance!` and $\mathbf{u_0}$ the initial storages (and the PID integrals which start out at $0$). +where $\mathbf{f}$ denotes `water_balance!` and $\mathbf{u_0} = \mathbf{0}$ the initial cumulative flows (and the PID integrals which also start out at $0$). -In general $\mathbf{f}$ is a non-linear function in $\mathbf{u}$. These non-linearities are introduced by: +In general $\mathbf{f}$ is a non-linear function in $\mathbf{u}$. These non-linearities are introduced by e.g.: - `ManningResistance` nodes; - `Basin` profiles; -- `TabulatedRatingCurve` Q(h) relations. +- `TabulatedRatingCurve` Q(h) relations +- `ContinuousControl` functions. The problem @eq-prob can be solved by various numerical time-integration methods. To do this the time interval $[t_0,t_\text{end}]$ is discretized into a finite number of time points $t_0 < t_1 < \ldots < t_N = t_\text{end}$ for which approximate solutions $\mathbf{w}_n \approx \mathbf{u}(t_n)$ are computed. In general we do not assume a fixed timestep (the interval between successive points in time). Rather, the solver attempts to make as large a step as possible while keeping error tolerances within requirements. The [solver settings](/reference/usage.qmd#sec-solver-settings) section details the available configuration options. diff --git a/python/ribasim_testmodels/ribasim_testmodels/backwater.py b/python/ribasim_testmodels/ribasim_testmodels/backwater.py index 04b79a550..b40b72083 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/backwater.py +++ b/python/ribasim_testmodels/ribasim_testmodels/backwater.py @@ -69,7 +69,7 @@ def backwater_model(): model.basin.add( Node(102, Point(1010.0, 0.0)), - [basin.State(level=[2.0]), basin.Profile(level=[0.0, 1.0], area=1e20)], + [basin.State(level=[2.0]), basin.Profile(level=[0.0, 1.0], area=1e10)], ) model.edge.add( model.manning_resistance[101],