Always require 'Basin / state' as an initial condition (#977)

Before we put a volume of 1 m3 in if no initial condition was given. But
since this will not always work well, and there is no other number that
makes a good default, it is best to remove this feature and request an
explicit initial condition.

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)
 

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
 
 """

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

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")

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.

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

python/ribasim_testmodels/ribasim_testmodels/backwater.py

@@ -52,24 +52,26 @@ 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,
             "precipitation": 0.0,
             "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(

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.

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.

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(

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(

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.