From 160b187da31f0c329e478820b304e7015680339f Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Jun 2023 15:30:31 +0200 Subject: [PATCH] More accurate interpolation --- core/src/bmi.jl | 3 +- core/src/create.jl | 3 +- core/src/io.jl | 3 +- core/src/solve.jl | 14 +-- core/src/utils.jl | 110 +++++++++++++++--- core/test/basin.jl | 45 ++++++- core/test/control.jl | 2 +- .../ribasim_testmodels/control.py | 2 +- 8 files changed, 148 insertions(+), 34 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index d3a32698e..57d0ddf83 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -144,7 +144,8 @@ function get_value(p::Parameters, node_id::Int, variable::String, u) # NOTE: Getting the level with get_level does NOT work since water_balance! # is not called during rootfinding for callback hasindex, basin_idx = id_index(basin.node_id, node_id) - value = basin.level[basin_idx](u[basin_idx]) + _, level = get_area_and_level(basin, basin_idx, u[basin_idx]) + value = level else throw(ValueError("Unsupported condition variable $variable.")) end diff --git a/core/src/create.jl b/core/src/create.jl index 5bb16ffe9..8a186f1e5 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -110,7 +110,7 @@ function Basin(db::DB, config::Config)::Basin infiltration = fill(NaN, length(node_id)) table = (; precipitation, potential_evaporation, drainage, infiltration) - area, level = create_storage_tables(db, config) + area, level, storage = create_storage_tables(db, config) # both static and forcing are optional, but we need fallback defaults static = load_structvector(db, config, BasinStaticV1) @@ -129,6 +129,7 @@ function Basin(db::DB, config::Config)::Basin current_level, area, level, + storage, time, ) end diff --git a/core/src/io.jl b/core/src/io.jl index 7b141f76a..e8ab71b90 100644 --- a/core/src/io.jl +++ b/core/src/io.jl @@ -131,7 +131,8 @@ function write_basin_output(model::Model) storage = reshape(vec(sol), nbasin, ntsteps) level = zero(storage) for (i, basin_storage) in enumerate(eachrow(storage)) - level[i, :] = p.basin.level[i].(basin_storage) + level[i, :] = + [get_area_and_level(p.basin, i, storage)[2] for storage in basin_storage] end basin = (; time, node_id, storage = vec(storage), level = vec(level)) diff --git a/core/src/solve.jl b/core/src/solve.jl index 7105dccb7..2191d2352 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -69,9 +69,10 @@ struct Basin{C} <: AbstractParameterNode infiltration::Vector{Float64} # cache this to avoid recomputation current_level::Vector{Float64} - # f(storage) - area::Vector{Interpolation} - level::Vector{Interpolation} + # Discrete values for interpolation + area::Vector{Vector{Float64}} + level::Vector{Vector{Float64}} + storage::Vector{Vector{Float64}} # data source for parameter updates time::StructVector{BasinForcingV1, C, Int} end @@ -243,11 +244,10 @@ Currently at less than 0.1 m. function formulate!(du::AbstractVector, basin::Basin, u::AbstractVector, t::Real)::Nothing for i in eachindex(du) storage = u[i] - area = basin.area[i](storage) - level = basin.level[i](storage) + area, level = get_area_and_level(basin, i, storage) basin.current_level[i] = level - bottom = basin_bottom_index(basin, i) - fixed_area = median(basin.area[i].u) + bottom = basin.level[i][1] + fixed_area = median(basin.area[i]) depth = max(level - bottom, 0.0) reduction_factor = min(depth, 0.1) / 0.1 diff --git a/core/src/utils.jl b/core/src/utils.jl index 255311c33..d702dd9a0 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -34,23 +34,103 @@ function profile_storage(levels::Vector{Float64}, areas::Vector{Float64})::Vecto return storages end -"Read the Basin / profile table and return all A(S) and h(S) functions" +"Read the Basin / profile table and return all area and level and computed storage values" function create_storage_tables( db::DB, config::Config, -)::Tuple{Vector{Interpolation}, Vector{Interpolation}} +)::Tuple{Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}} profiles = load_structvector(db, config, BasinProfileV1) - area = Interpolation[] - level = Interpolation[] + area = Vector{Vector{Float64}}() + level = Vector{Vector{Float64}}() + storage = Vector{Vector{Float64}}() + for group in IterTools.groupby(row -> row.node_id, profiles) group_area = getproperty.(group, :area) group_level = getproperty.(group, :level) group_storage = profile_storage(group_level, group_area) - area_itp = LinearInterpolation(group_area, group_storage) - level_itp = LinearInterpolation(group_level, group_storage) - push!(area, area_itp) - push!(level, level_itp) + push!(area, group_area) + push!(level, group_level) + push!(storage, group_storage) + end + return area, level, storage +end + +function get_area_and_level( + basin::Basin, + state_idx::Int, + storage::Float64, +)::Tuple{Float64, Float64} + storage_discrete = basin.storage[state_idx] + area_discrete = basin.area[state_idx] + level_discrete = basin.level[state_idx] + + return get_area_and_level(storage_discrete, area_discrete, level_discrete, storage) +end + +function get_area_and_level( + storage_discrete::Vector{Float64}, + area_discrete::Vector{Float64}, + level_discrete::Vector{Float64}, + storage::Float64, +)::Tuple{Float64, Float64} + # storage_idx: smallest index such that storage_discrete[storage_idx] >= storage + storage_idx = searchsortedfirst(storage_discrete, storage) + + if storage_idx == 1 # If the lowest discrete_storage level 0, this can only happen if the storage is 0 since storage is never negative + level = level_discrete[1] + area = area_discrete[1] + + elseif storage_idx == length(storage_discrete) + 1 + # These values yield linear extrapolation of area(level) based on the last 2 values + area_lower = area_discrete[end - 1] + area_higher = area_discrete[end] + level_lower = level_discrete[end - 1] + level_higher = level_discrete[end] + storage_lower = storage_discrete[end - 1] + storage_higher = storage_discrete[end] + + area_diff = area_higher - area_lower + level_diff = level_higher - level_lower + + if area_diff ≈ 0 + # Constant area means linear interpolation of level + area = area_lower + + level = + area_higher + + level_diff * (storage - storage_higher) / (storage_higher - storage_lower) + else + area = sqrt( + area_higher^2 + 2 * (storage - storage_higher) * area_diff / level_diff, + ) + level = level_lower + level_diff * (area - area_lower) / area_diff + end + + else + area_lower = area_discrete[storage_idx - 1] + area_higher = area_discrete[storage_idx] + level_lower = level_discrete[storage_idx - 1] + level_higher = level_discrete[storage_idx] + storage_lower = storage_discrete[storage_idx - 1] + storage_higher = storage_discrete[storage_idx] + + area_diff = area_higher - area_lower + level_diff = level_higher - level_lower + + if area_diff ≈ 0 + # Constant area means linear interpolation of level + area = area_lower + level = + level_lower + + level_diff * (storage - storage_lower) / (storage_higher - storage_lower) + + else + area = + sqrt(area_lower^2 + 2 * (storage - storage_lower) * area_diff / level_diff) + level = level_lower + level_diff * (area - area_lower) / area_diff + end end + return area, level end @@ -206,23 +286,15 @@ function id_index(ids::Indices{Int}, id::Int) return hasindex, idx end -"Return the bottom elevation of the basin with index i" -function basin_bottom_index(basin::Basin, i::Int)::Float64 - # get level(storage) interpolation function - itp = basin.level[i] - # and return the first level in the underlying table, which represents the bottom - return first(itp.u) -end - "Return the bottom elevation of the basin with index i, or nothing if it doesn't exist" function basin_bottom(basin::Basin, node_id::Int)::Union{Float64, Nothing} basin = Dictionary(basin.node_id, basin.level) hasindex, token = gettoken(basin, node_id) return if hasindex # get level(storage) interpolation function - itp = gettokenvalue(basin, token) - # and return the first level in the underlying table, which represents the bottom - first(itp.u) + level_discrete = gettokenvalue(basin, token) + # and return the first level in this vector, representing the bottom + first(level_discrete) else nothing end diff --git a/core/test/basin.jl b/core/test/basin.jl index ca2494d11..f90a99530 100644 --- a/core/test/basin.jl +++ b/core/test/basin.jl @@ -17,7 +17,7 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test model.integrator.sol.retcode == Ribasim.ReturnCode.Success - @test model.integrator.sol.u[end] ≈ Float32[476.8749, 476.85898, 1.8706497, 786.96686] skip = + @test model.integrator.sol.u[end] ≈ Float32[452.9688, 453.0431, 1.8501105, 1238.0144] skip = Sys.isapple() end @@ -28,7 +28,7 @@ end @test model isa Ribasim.Model @test model.integrator.sol.retcode == Ribasim.ReturnCode.Success @test length(model.integrator.p.basin.precipitation) == 4 - @test model.integrator.sol.u[end] ≈ Float32[480.5936, 480.58762, 1.4178488, 793.6097] skip = + @test model.integrator.sol.u[end] ≈ Float32[428.06897, 428.07315, 1.3662858, 1249.2343] skip = Sys.isapple() end @@ -39,7 +39,46 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test model.integrator.sol.retcode == Ribasim.ReturnCode.Success - @test model.integrator.sol.u[end] ≈ Float32[27.273159, 340.68964] skip = Sys.isapple() + @test model.integrator.sol.u[end] ≈ Float32[1.4875988, 366.4851] 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.tables[end].t[end] == 1.2 end + +@testset "Profile" begin + n_interpolations = 100 + storage = range(0.0, 1000.0, n_interpolations) + + # Covers interpolation for constant and non-constant area, extrapolation for constant area + area_discrete = [0.0, 100.0, 100.0] + level_discrete = [0.0, 10.0, 15.0] + storage_discrete = Ribasim.profile_storage(level_discrete, area_discrete) + + area, level = zip( + [ + Ribasim.get_area_and_level(storage_discrete, area_discrete, level_discrete, s) for s in storage + ]..., + ) + + level_expected = + ifelse.(storage .< 500.0, sqrt.(storage ./ 5), 10.0 .+ (storage .- 500.0) ./ 100.0) + + @test all(level .≈ level_expected) + area_expected = min.(10.0 * level_expected, 100.0) + @test all(area .≈ area_expected) + + # Covers extrapolation for non-constant area + area_discrete = [0.0, 100.0] + level_discrete = [0.0, 10.0] + storage_discrete = Ribasim.profile_storage(level_discrete, area_discrete) + + area, level = zip( + [ + Ribasim.get_area_and_level(storage_discrete, area_discrete, level_discrete, s) for s in storage + ]..., + ) + + level_expected = sqrt.(storage ./ 5) + @test all(level .≈ level_expected) + area_expected = 10.0 * level_expected + @test all(area .≈ area_expected) +end diff --git a/core/test/control.jl b/core/test/control.jl index 3401a9263..b37e8a80b 100644 --- a/core/test/control.jl +++ b/core/test/control.jl @@ -1,6 +1,6 @@ import Ribasim -toml_path = normpath(@__DIR__, "../../data/pump_control/control.toml") +toml_path = normpath(@__DIR__, "../../data/pump_control/pump_control.toml") @testset "pump control" begin @test ispath(toml_path) diff --git a/python/ribasim_testmodels/ribasim_testmodels/control.py b/python/ribasim_testmodels/ribasim_testmodels/control.py index 94cd71f8f..374fea6fe 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/control.py @@ -129,7 +129,7 @@ def pump_control_model() -> ribasim.Model: # Setup a model: model = ribasim.Model( - modelname="control", + modelname="pump_control", node=node, edge=edge, basin=basin,