From 75635c9b71fd90b62633ce21572625ee630a8ea1 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 22 Jan 2024 11:17:20 +0100 Subject: [PATCH 1/2] Always use initial conditions --- core/test/bmi_test.jl | 2 +- core/test/run_models_test.jl | 2 +- .../ribasim_testmodels/backwater.py | 8 +++-- .../ribasim_testmodels/basic.py | 11 +++++-- .../ribasim_testmodels/discrete_control.py | 4 ++- .../ribasim_testmodels/invalid.py | 31 ++++++++++++++++--- .../ribasim_testmodels/time.py | 9 +++++- .../ribasim_testmodels/trivial.py | 4 ++- 8 files changed, 56 insertions(+), 15 deletions(-) diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index 8c6c95dc7..39599eeac 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -49,7 +49,7 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) storage0 = BMI.get_value_ptr(model, "volume") - @test storage0 == ones(4) + @test storage0 ≈ ones(4) @test_throws "Unknown variable foo" BMI.get_value_ptr(model, "foo") BMI.update_until(model, 86400.0) storage = BMI.get_value_ptr(model, "volume") diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 995b07032..227ec010f 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -73,7 +73,7 @@ @test flow.from_node_id[1:4] == [6, typemax(Int), 0, 6] @test flow.to_node_id[1:4] == [6, typemax(Int), typemax(Int), 0] - @test basin.storage[1] == 1.0 + @test basin.storage[1] ≈ 1.0 @test basin.level[1] ≈ 0.044711584 # The exporter interpolates 1:1 for three subgrid elements, but shifted by 1.0 meter. diff --git a/python/ribasim_testmodels/ribasim_testmodels/backwater.py b/python/ribasim_testmodels/ribasim_testmodels/backwater.py index e9adc4f07..42e950ec9 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/backwater.py +++ b/python/ribasim_testmodels/ribasim_testmodels/backwater.py @@ -52,16 +52,17 @@ def backwater_model(): ) # Rectangular profile, width of 1.0 m. + basin_ids = ids[node_type == "Basin"] profile = pd.DataFrame( data={ - "node_id": np.repeat(ids[node_type == "Basin"], 2), + "node_id": np.repeat(basin_ids, 2), "area": [20.0, 20.0] * n_basin, "level": [0.0, 1.0] * n_basin, } ) static = pd.DataFrame( data={ - "node_id": ids[node_type == "Basin"], + "node_id": basin_ids, "drainage": 0.0, "potential_evaporation": 0.0, "infiltration": 0.0, @@ -69,7 +70,8 @@ def backwater_model(): "urban_runoff": 0.0, } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame(data={"node_id": basin_ids, "level": 0.05}) + basin = ribasim.Basin(profile=profile, static=static, state=state) manning_resistance = ribasim.ManningResistance( static=pd.DataFrame( diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 7a3982173..657314b95 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -36,8 +36,10 @@ def basic_model() -> ribasim.Model: ) static = static.iloc[[0, 0, 0, 0]] static["node_id"] = [1, 3, 6, 9] - - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame( + data={"node_id": static["node_id"], "level": 0.04471158417652035} + ) + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup linear resistance: linear_resistance = ribasim.LinearResistance( @@ -303,8 +305,11 @@ def tabulated_rating_curve_model() -> ribasim.Model: "urban_runoff": 0.0, } ) + state = pd.DataFrame( + data={"node_id": static["node_id"], "level": 0.04471158417652035} + ) - basin = ribasim.Basin(profile=profile, static=static) + basin = ribasim.Basin(profile=profile, static=static, state=state) # Set up a rating curve node: # Discharge: lose 1% of storage volume per day at storage = 1000.0. diff --git a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py index 6481214ad..e82b606a0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py @@ -511,7 +511,9 @@ def tabulated_rating_curve_control_model() -> ribasim.Model: } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame(data={"node_id": [1], "level": 0.04471158417652035}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) # Set up a rating curve node: # Discharge: lose 1% of storage volume per day at storage = 100.0. diff --git a/python/ribasim_testmodels/ribasim_testmodels/invalid.py b/python/ribasim_testmodels/ribasim_testmodels/invalid.py index 19cee6613..7301f7f2a 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/invalid.py +++ b/python/ribasim_testmodels/ribasim_testmodels/invalid.py @@ -61,7 +61,9 @@ def invalid_qh_model(): } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame(data={"node_id": [3], "level": 1.4112729908597084}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) rating_curve_static = pd.DataFrame( # Invalid: levels must not be repeated @@ -171,7 +173,14 @@ def invalid_fractional_flow_model(): } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame( + data={ + "node_id": [1, 2], + "level": 1.4112729908597084, + } + ) + + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup terminal: terminal = ribasim.Terminal(static=pd.DataFrame(data={"node_id": [5, 6]})) @@ -263,7 +272,14 @@ def invalid_discrete_control_model(): } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame( + data={ + "node_id": [1, 3], + "level": 1.4112729908597084, + } + ) + + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup pump: pump = ribasim.Pump( @@ -391,7 +407,14 @@ def invalid_edge_types_model(): } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame( + data={ + "node_id": [1, 3], + "level": 0.04471158417652035, + } + ) + + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup pump: pump = ribasim.Pump( diff --git a/python/ribasim_testmodels/ribasim_testmodels/time.py b/python/ribasim_testmodels/ribasim_testmodels/time.py index bb34dfea7..42a0bdd6a 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/time.py +++ b/python/ribasim_testmodels/ribasim_testmodels/time.py @@ -66,7 +66,14 @@ def flow_boundary_time_model(): } ) - basin = ribasim.Basin(profile=profile, static=static) + state = pd.DataFrame( + data={ + "node_id": [2], + "level": 0.04471158417652035, + } + ) + + basin = ribasim.Basin(profile=profile, static=static, state=state) n_times = 100 time = pd.date_range( diff --git a/python/ribasim_testmodels/ribasim_testmodels/trivial.py b/python/ribasim_testmodels/ribasim_testmodels/trivial.py index 664b69197..99b5067cf 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/trivial.py +++ b/python/ribasim_testmodels/ribasim_testmodels/trivial.py @@ -74,6 +74,8 @@ def trivial_model() -> ribasim.Model: } ) + state = pd.DataFrame(data={"node_id": [6], "level": 0.04471158417652035}) + # Create a subgrid level interpolation from one basin to three elements. Scale one to one, but: # # 22. start at -1.0 @@ -88,7 +90,7 @@ def trivial_model() -> ribasim.Model: "subgrid_level": [-1.0, 0.0, 0.0, 1.0, 1.0, 2.0], } ) - basin = ribasim.Basin(profile=profile, static=static, subgrid=subgrid) + basin = ribasim.Basin(profile=profile, static=static, state=state, subgrid=subgrid) # Set up a rating curve node: # Discharge: lose 1% of storage volume per day at storage = 1000.0. From 6dee56358f355dc4f8aff24fa9ca7a03d7d3059a Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 22 Jan 2024 11:39:02 +0100 Subject: [PATCH 2/2] finish --- core/src/bmi.jl | 12 ++---------- core/src/utils.jl | 26 +++++++++++++++++++------- core/src/validation.jl | 24 ------------------------ core/test/utils_test.jl | 4 +--- 4 files changed, 22 insertions(+), 44 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index d2d7f43bb..7f11ee34d 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -73,16 +73,8 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model end @debug "Read database into memory." - storage = if isempty(state) - # default to nearly empty basins, perhaps make required input - fill(1.0, n) - else - storages, errors = get_storages_from_levels(parameters.basin, state.level) - if errors - error("Encountered errors while parsing the initial levels of basins.") - end - storages - end + storage = get_storages_from_levels(parameters.basin, state.level) + # Synchronize level with storage set_current_basin_properties!(parameters.basin, storage) diff --git a/core/src/utils.jl b/core/src/utils.jl index 5040f786c..20ff9f75c 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -374,16 +374,28 @@ function get_storage_from_level(basin::Basin, state_idx::Int, level::Float64)::F end """Compute the storages of the basins based on the water level of the basins.""" -function get_storages_from_levels( - basin::Basin, - levels::Vector, -)::Tuple{Vector{Float64}, Bool} - storages = Float64[] +function get_storages_from_levels(basin::Basin, levels::Vector)::Vector{Float64} + errors = false + state_length = length(levels) + basin_length = length(basin.level) + if state_length != basin_length + @error "Unexpected 'Basin / state' length." state_length basin_length + errors = true + end + storages = zeros(state_length) for (i, level) in enumerate(levels) - push!(storages, get_storage_from_level(basin, i, level)) + storage = get_storage_from_level(basin, i, level) + if isnan(storage) + errors = true + end + storages[i] = storage + end + if errors + error("Encountered errors while parsing the initial levels of basins.") end - return storages, any(isnan.(storages)) + + return storages end """ diff --git a/core/src/validation.jl b/core/src/validation.jl index d49dc3506..922b02fe2 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -338,30 +338,6 @@ function variable_nt(s::Any) NamedTuple{names}((getfield(s, x) for x in names)) end -function is_consistent(node, edge, state, static, profile, time) - - # Check that node ids exist - # TODO Do we need to check the reverse as well? All ids in use? - ids = node.fid - @assert edge.from_node_id ⊆ ids "Edge from_node_id not in node ids" - @assert edge.to_node_id ⊆ ids "Edge to_node_id not in node ids" - @assert state.node_id ⊆ ids "State id not in node ids" - @assert static.node_id ⊆ ids "Static id not in node ids" - @assert profile.node_id ⊆ ids "Profile id not in node ids" - @assert time.node_id ⊆ ids "Time id not in node ids" - - # Check edges for uniqueness - @assert allunique(edge, [:from_node_id, :to_node_id]) "Duplicate edge found" - - # TODO Check states - - # TODO Check statics - - # TODO Check forcings - - true -end - # functions used by sort(x; by) sort_by_fid(row) = row.fid sort_by_id(row) = row.node_id diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index ed2bc1a1e..d4263a878 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -125,9 +125,7 @@ end logger = TestLogger() with_logger(logger) do - storages, errors = Ribasim.get_storages_from_levels(basin, [-1.0]) - @test isnan(storages[1]) - @test errors + @test_throws ErrorException Ribasim.get_storages_from_levels(basin, [-1.0]) end @test length(logger.logs) == 1