From fb54b003a38e0af71978857ae77c24fb75695dcb Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Wed, 14 Jun 2023 11:23:50 +0200 Subject: [PATCH] calculate storage from Basin / profile (#280) Fixes #225 This removes the need to provide the storage for a profile. We can calculate it ourselves by integrating over A(h) to get A(S) and h(S). First this updates the test models to a correct storage in the profile. Then we calculate the profile and get the same results. Then we remove the storage from the schemas and docs. I hoped this would help with #80 but based on the plot I added there it doesn't seem to. This does not yet add all the validation and extrapolation for lookup tables as described in #279, that can be done separately. --- .github/pull_request_template.md | 2 +- core/src/Ribasim.jl | 2 +- core/src/config.jl | 23 ++++++++++ core/src/create.jl | 16 ------- core/src/utils.jl | 46 ++++++++++++------- core/src/validation.jl | 9 ++-- core/test/basin.jl | 6 +-- core/test/utils.jl | 7 +++ docs/contribute/addnode.qmd | 2 +- docs/core/usage.qmd | 22 +++++---- docs/python/examples.ipynb | 4 +- docs/schema/BasinProfile.schema.json | 6 --- python/ribasim/ribasim/models.py | 1 - python/ribasim/ribasim/node_types/basin.py | 4 +- .../ribasim_testmodels/basic.py | 4 +- qgis/core/nodes.py | 1 - 16 files changed, 87 insertions(+), 68 deletions(-) diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 4ec7ee803..e2e198a22 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,4 +1,4 @@ -Solves issue #? +Fixes # # Description diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 564b1968d..4fdafb762 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -27,8 +27,8 @@ TimerOutputs.complement!() include("validation.jl") include("solve.jl") -include("utils.jl") include("config.jl") +include("utils.jl") include("lib.jl") include("io.jl") include("create.jl") diff --git a/core/src/config.jl b/core/src/config.jl index 0e63986f1..7e3083125 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -12,6 +12,29 @@ for sv in nodeschemas push!(nodekinds[node], kind) end +""" +Add fieldnames with Maybe{String} type to struct expression. Requires @option use before it. +""" +macro addfields(typ::Expr, fieldnames) + for fieldname in fieldnames + push!(typ.args[3].args, Expr(:(=), Expr(:(::), fieldname, Maybe{String}), nothing)) + end + return esc(typ) +end + +""" +Add all TableOption subtypes as fields to struct expression. Requires @option use before it. +""" +macro addnodetypes(typ::Expr) + for nodetype in nodetypes + push!( + typ.args[3].args, + Expr(:(=), Expr(:(::), nodetype, nodetype), Expr(:call, nodetype)), + ) + end + return esc(typ) +end + # Generate structs for each nodetype for use in Config abstract type TableOption end for (T, kinds) in pairs(nodekinds) diff --git a/core/src/create.jl b/core/src/create.jl index a1445e1f5..5bb16ffe9 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -52,22 +52,6 @@ function ManningResistance(db::DB, config::Config)::ManningResistance ) end -function create_storage_tables(db::DB, config::Config) - profiles = load_structvector(db, config, BasinProfileV1) - area = Interpolation[] - level = Interpolation[] - for group in IterTools.groupby(row -> row.node_id, profiles) - group_storage = getproperty.(group, :storage) - group_area = getproperty.(group, :area) - group_level = getproperty.(group, :level) - area_itp = LinearInterpolation(group_area, group_storage) - level_itp = LinearInterpolation(group_level, group_storage) - push!(area, area_itp) - push!(level, level_itp) - end - return area, level -end - function FractionalFlow(db::DB, config::Config)::FractionalFlow static = load_structvector(db, config, FractionalFlowStaticV1) return FractionalFlow(static.node_id, static.fraction) diff --git a/core/src/utils.jl b/core/src/utils.jl index 9814e7688..9dc9f149e 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -19,27 +19,39 @@ function create_graph( return graph, edge_ids, edge_connection_types end -""" -Add fieldnames with Maybe{String} type to struct expression. Requires @option use before it. -""" -macro addfields(typ::Expr, fieldnames) - for fieldname in fieldnames - push!(typ.args[3].args, Expr(:(=), Expr(:(::), fieldname, Maybe{String}), nothing)) +"Calculate a profile storage by integrating the areas over the levels" +function profile_storage(levels::Vector{Float64}, areas::Vector{Float64})::Vector{Float64} + # profile starts at the bottom; first storage is 0 + storages = zero(areas) + n = length(storages) + + for i in 2:n + Δh = levels[i] - levels[i - 1] + avg_area = 0.5 * (areas[i - 1] + areas[i]) + ΔS = avg_area * Δh + storages[i] = storages[i - 1] + ΔS end - return esc(typ) + return storages end -""" -Add all TableOption subtypes as fields to struct expression. Requires @option use before it. -""" -macro addnodetypes(typ::Expr) - for nodetype in nodetypes - push!( - typ.args[3].args, - Expr(:(=), Expr(:(::), nodetype, nodetype), Expr(:call, nodetype)), - ) +"Read the Basin / profile table and return all A(S) and h(S) functions" +function create_storage_tables( + db::DB, + config::Config, +)::Tuple{Vector{Interpolation}, Vector{Interpolation}} + profiles = load_structvector(db, config, BasinProfileV1) + area = Interpolation[] + level = Interpolation[] + 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) end - return esc(typ) + return area, level end """ diff --git a/core/src/validation.jl b/core/src/validation.jl index 4b31461f1..f62551383 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -84,7 +84,6 @@ end @version BasinProfileV1 begin node_id::Int - storage::Float64 area::Float64 level::Float64 end @@ -195,19 +194,21 @@ end sort_by_id(row) = row.node_id sort_by_time_id(row) = (row.time, row.node_id) sort_by_id_level(row) = (row.node_id, row.level) -sort_by_id_storage(row) = (row.node_id, row.storage) # get the right sort by function given the Schema, with sort_by_id as the default sort_by_function(table::StructVector{<:Legolas.AbstractRecord}) = sort_by_id -sort_by_function(table::StructVector{TabulatedRatingCurveStaticV1}) = sort_by_id_level -sort_by_function(table::StructVector{BasinProfileV1}) = sort_by_id_storage const TimeSchemas = Union{TabulatedRatingCurveTimeV1, BasinForcingV1} +const LevelLookupSchemas = Union{TabulatedRatingCurveStaticV1, BasinProfileV1} function sort_by_function(table::StructVector{<:TimeSchemas}) return sort_by_time_id end +function sort_by_function(table::StructVector{<:LevelLookupSchemas}) + return sort_by_id_level +end + """ Depending on if a table can be sorted, either sort it or assert that it is sorted. diff --git a/core/test/basin.jl b/core/test/basin.jl index 8df74e1ad..f399dc432 100644 --- a/core/test/basin.jl +++ b/core/test/basin.jl @@ -9,7 +9,7 @@ using SciMLBase 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[948.93713, 948.96857, 1.8637465, 1573.5077] skip = + @test model.integrator.sol.u[end] ≈ Float32[476.8749, 476.85898, 1.8706497, 786.96686] skip = Sys.isapple() end @@ -20,7 +20,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[915.72656, 915.701, 1.3522862, 1583.116] skip = + @test model.integrator.sol.u[end] ≈ Float32[480.5936, 480.58762, 1.4178488, 793.6097] skip = Sys.isapple() end @@ -31,7 +31,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[54.459435, 313.50992] skip = Sys.isapple() + @test model.integrator.sol.u[end] ≈ Float32[27.273159, 340.68964] 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 diff --git a/core/test/utils.jl b/core/test/utils.jl index 9c55d925c..f3bf5fdec 100644 --- a/core/test/utils.jl +++ b/core/test/utils.jl @@ -7,3 +7,10 @@ using Test @test Ribasim.id_index(ids, 4) === (true, 2) @test Ribasim.id_index(ids, 5) === (false, 0) end + +@testset "profile_storage" begin + @test Ribasim.profile_storage([0.0, 1.0], [0.0, 1000.0]) == [0.0, 500.0] + @test Ribasim.profile_storage([6.0, 7.0], [0.0, 1000.0]) == [0.0, 500.0] + @test Ribasim.profile_storage([6.0, 7.0, 9.0], [0.0, 1000.0, 1000.0]) == + [0.0, 500.0, 2500.0] +end diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd index bf753195d..5ec0d26c9 100644 --- a/docs/contribute/addnode.qmd +++ b/docs/contribute/addnode.qmd @@ -184,7 +184,7 @@ Future: Test for the specific behaviour of each node type There are several parts of the documentation which should be updated with the new node type: - `docs/core/equations` should contain a short explanation and if possible an analytical expression for the behaviour of the new node; -- `docs/core/usage.qmd` should conntain a short explanation of the node and the possible schemas associated with it; +- `docs/core/usage.qmd` should contain a short explanation of the node and the possible schemas associated with it; - The example models constructed in `docs/python/examples.ipynb` should be extended with the new node type or a new example model with the new node type should be made. # Finishing up diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 93e302ec0..267b8e762 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -235,22 +235,24 @@ The profile table defines the physical dimensions of the storage reservoir of ea column | type | unit | restriction --------- | ------- | ------------ | ----------- node_id | Int | - | sorted -storage | Float64 | $m^3$ | non-negative, per node_id: start at 0 and increasing -area | Float64 | $m^2$ | non-negative -level | Float64 | $m$ | - +area | Float64 | $m^2$ | non-negative, per node_id: start at 0 and increasing +level | Float64 | $m$ | per node_id: increasing The level is the level at the basin outlet. All levels are defined in meters above a datum that is the same for the entire model. An example of the first 5 rows of such a table is given below. The first 4 rows define the profile of ID `2`. The number of rows can vary per ID. Using a very large number of rows may impact performance. -node_id | storage | area | level -------- |--------------- |-------------- |------- - 2 | 0.0 | 1.36404e5 | -0.105 - 2 | 24726.2 | 1.36404e5 | 0.095 - 2 | 49452.5 | 1.36404e5 | 0.295 - 2 | 2.49735e6 | 1.36404e5 | 20.095 - 3 | 0.0 | 50663.3 | 2.129 +node_id | area | level +------- |------- |------- + 2 | 0.0 | 6.0 + 2 | 1000.0 | 7.0 + 2 | 1000.0 | 9.0 + 3 | 0.0 | 2.2 + +We use the symbol $A$ for area, $h$ for level and $S$ for storage. +The profile provides a function $A(h)$ for each basin. +Internally this get converted to two functions, $A(S)$ and $h(S)$, by integrating over the function, setting the storage to zero for the bottom of the profile. ### FractionalFlow diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index b24dd03ae..5f61d7506 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -162,7 +162,7 @@ "profile = pd.DataFrame(\n", " data={\n", " \"node_id\": [1, 1, 3, 3, 6, 6, 9, 9],\n", - " \"storage\": [0.0, 1000.0] * 4,\n", + " \"storage\": [0.0, 500.0] * 4,\n", " \"area\": [0.0, 1000.0] * 4,\n", " \"level\": [0.0, 1.0] * 4,\n", " }\n", @@ -1000,7 +1000,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.10.11" } }, "nbformat": 4, diff --git a/docs/schema/BasinProfile.schema.json b/docs/schema/BasinProfile.schema.json index d3ba6fa22..9f771318f 100644 --- a/docs/schema/BasinProfile.schema.json +++ b/docs/schema/BasinProfile.schema.json @@ -12,11 +12,6 @@ "description": "area", "type": "number" }, - "storage": { - "format": "double", - "description": "storage", - "type": "number" - }, "node_id": { "format": "default", "description": "node_id", @@ -30,7 +25,6 @@ }, "required": [ "node_id", - "storage", "area", "level" ], diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index 1cd1b8d97..873fd2ee5 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -119,7 +119,6 @@ class BasinState(BaseModel): class BasinProfile(BaseModel): remarks: Optional[str] = Field("", description="a hack for pandera") area: float = Field(..., description="area") - storage: float = Field(..., description="storage") node_id: int = Field(..., description="node_id") level: float = Field(..., description="level") diff --git a/python/ribasim/ribasim/node_types/basin.py b/python/ribasim/ribasim/node_types/basin.py index 1d3460905..147524532 100644 --- a/python/ribasim/ribasim/node_types/basin.py +++ b/python/ribasim/ribasim/node_types/basin.py @@ -64,9 +64,7 @@ class Config: validate_assignment = True def sort(self): - self.profile = self.profile.sort_values( - ["node_id", "storage"], ignore_index=True - ) + self.profile = self.profile.sort_values(["node_id", "level"], ignore_index=True) if self.static is not None: self.static = self.static.sort_values("node_id", ignore_index=True) if self.forcing is not None: diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 22a58ae68..b19ac8204 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -85,7 +85,7 @@ def basic_model() -> ribasim.Model: profile = pd.DataFrame( data={ "node_id": [1, 1, 3, 3, 6, 6, 9, 9], - "storage": [0.0, 1000.0] * 4, + "storage": [0.0, 500.0] * 4, "area": [0.0, 1000.0] * 4, "level": [0.0, 1.0] * 4, } @@ -323,7 +323,7 @@ def tabulated_rating_curve_model() -> ribasim.Model: profile = pd.DataFrame( data={ "node_id": [1, 1, 4, 4], - "storage": [0.0, 1000.0] * 2, + "storage": [0.0, 500.0] * 2, "area": [0.0, 1000.0] * 2, "level": [0.0, 1.0] * 2, } diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py index 0d3e80363..f5fadffae 100644 --- a/qgis/core/nodes.py +++ b/qgis/core/nodes.py @@ -262,7 +262,6 @@ class BasinProfile(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), - QgsField("storage", QVariant.Double), QgsField("area", QVariant.Double), QgsField("level", QVariant.Double), ]