require profile areas to start positive (#475)

Fixes #472 Since we change the validation from requiring profile areas to start at 0 to requiring a positive area, this breaks all existing models. So therefore all testmodels also had to be updated. I changed 0 to 0.01 m2 (10 x 10cm), sufficiently small not to have to update test states. Or for the cases where we added a sharply widening profile from level 0 to 1 followed by a vertical profile, I just make the profile start at the intended value. `get_area_and_level` did not need updating since 0 was not hardcoded. There is one unrelated change with regard to the logging of a vector of error messages, I change that since I think it now looks clearer: ```julia julia> errors = ["this is wrong", "that is wrong"] julia> @error join(errors, "\n") ┌ Error: this is wrong │ that is wrong └ @ Main REPL[2]:1 julia> foreach(x -> @error(x), errors) ┌ Error: this is wrong └ @ Main REPL[3]:1 ┌ Error: that is wrong └ @ Main REPL[3]:1 ```
Deltares · Aug 3, 2023 · b87d90f · b87d90f
1 parent daeb7bd
commit b87d90f
Show file tree

Hide file tree

Showing 15 changed files with 41 additions and 48 deletions.
diff --git a/core/src/solve.jl b/core/src/solve.jl
@@ -118,7 +118,7 @@ struct Basin{C} <: AbstractParameterNode
                 time,
             )
         else
-            @error join(errors, "\n")
+            foreach(x -> @error(x), errors)
             error("Errors occurred when parsing Basin data.")
         end
     end
@@ -349,7 +349,7 @@ function valid_n_neighbors(p::Parameters)::Bool
     if isempty(errors)
         return true
     else
-        @error join(errors, "\n")
+        foreach(x -> @error(x), errors)
         return false
     end
 end

diff --git a/core/src/validation.jl b/core/src/validation.jl
@@ -255,8 +255,6 @@ function is_consistent(node, edge, state, static, profile, forcing)
 
     # TODO Check statics
 
-    # TODO Check profiles
-
     # TODO Check forcings
 
     true
@@ -323,13 +321,13 @@ function valid_edges(
     if isempty(errors)
         return true
     else
-        @error join(errors, "\n")
+        foreach(x -> @error(x), errors)
         return false
     end
 end
 
 """
-Check whether the profile data has no repeats in the levels and the areas start at 0.
+Check whether the profile data has no repeats in the levels and the areas start positive.
 """
 function valid_profiles(
     node_id::Indices{Int},
@@ -343,10 +341,10 @@ function valid_profiles(
             push!(errors, "Basin #$id has repeated levels, this cannot be interpolated.")
         end
 
-        if areas[1] != 0
+        if areas[1] <= 0
             push!(
                 errors,
-                "Basin profiles must start with area 0 at the bottom (got area $(areas[1]) for node #$id).",
+                "Basin profiles cannot start with area <= 0 at the bottom for numerical reasons (got area $(areas[1]) for node #$id).",
             )
         end
     end

diff --git a/core/test/utils.jl b/core/test/utils.jl
@@ -20,7 +20,7 @@ end
 
 @testset "bottom" begin
     # create two basins with different bottoms/levels
-    area = [[0.0, 1.0], [0.0, 1.0]]
+    area = [[0.01, 1.0], [0.01, 1.0]]
     level = [[0.0, 1.0], [4.0, 5.0]]
     storage = Ribasim.profile_storage.(level, area)
     target_level = [0.0, 0.0]

diff --git a/core/test/validation.jl b/core/test/validation.jl
@@ -8,10 +8,10 @@ using Logging
 @testset "Basin profile validation" begin
     node_id = Indices([1])
     level = [[0.0, 0.0]]
-    area = [[100.0, 100.0]]
+    area = [[0.0, 100.0]]
     errors = Ribasim.valid_profiles(node_id, level, area)
     @test "Basin #1 has repeated levels, this cannot be interpolated." in errors
-    @test "Basin profiles must start with area 0 at the bottom (got area 100.0 for node #1)." in
+    @test "Basin profiles cannot start with area <= 0 at the bottom for numerical reasons (got area 0.0 for node #1)." in
           errors
     @test length(errors) == 2
 

diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd
@@ -207,7 +207,7 @@ The profile table defines the physical dimensions of the storage reservoir of ea
 column    | type    | unit         | restriction
 --------- | ------- | ------------ | -----------
 node_id   | Int     | -            | sorted
-area      | Float64 | $m^2$        | non-negative, per node_id: start at 0 and increasing
+area      | Float64 | $m^2$        | non-negative, per node_id: start positive 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
@@ -217,14 +217,15 @@ per ID. Using a very large number of rows may impact performance.
 
 node_id |   area | level
 ------- |------- |-------
-      2 |    0.0 |   6.0
+      2 |    1.0 |   6.0
       2 | 1000.0 |   7.0
       2 | 1000.0 |   9.0
-      3 |    0.0 |   2.2
+      3 |    1.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.
+The minimum area cannot be zero to avoid numerical issues.
 The maximum area is used to convert the precipitation flux into an inflow.
 
 ### FractionalFlow

diff --git 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",
-    "        \"area\": [0.0, 1000.0] * 4,\n",
+    "        \"area\": [0.01, 1000.0] * 4,\n",
     "        \"level\": [0.0, 1.0] * 4,\n",
     "    }\n",
     ")\n",
@@ -736,9 +736,9 @@
    "source": [
     "profile = pd.DataFrame(\n",
     "    data={\n",
-    "        \"node_id\": [1, 1, 1, 3, 3, 3],\n",
-    "        \"area\": [0.0, 100.0, 100.0] * 2,\n",
-    "        \"level\": [0.0, 0.001, 1.0] * 2,\n",
+    "        \"node_id\": [1, 1, 3, 3],\n",
+    "        \"area\": [100.0, 100.0] * 2,\n",
+    "        \"level\": [0.0, 1.0] * 2,\n",
     "    }\n",
     ")\n",
     "\n",
@@ -980,13 +980,6 @@
    "source": [
     "model.print_discrete_control_record(datadir / \"control/output/control.arrow\")"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

diff --git a/python/ribasim_testmodels/ribasim_testmodels/backwater.py b/python/ribasim_testmodels/ribasim_testmodels/backwater.py
@@ -54,9 +54,9 @@ def backwater_model():
     # Rectangular profile, width of 1.0 m.
     profile = pd.DataFrame(
         data={
-            "node_id": np.repeat(ids[node_type == "Basin"], 3),
-            "area": [0.0, 20.0, 20.0] * n_basin,
-            "level": [0.0, 0.01, 1.0] * n_basin,
+            "node_id": np.repeat(ids[node_type == "Basin"], 2),
+            "area": [20.0, 20.0] * n_basin,
+            "level": [0.0, 1.0] * n_basin,
         }
     )
     static = pd.DataFrame(

diff --git 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],
-            "area": [0.0, 1000.0] * 4,
+            "area": [0.01, 1000.0] * 4,
             "level": [0.0, 1.0] * 4,
         }
     )
@@ -322,7 +322,7 @@ def tabulated_rating_curve_model() -> ribasim.Model:
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1, 4, 4],
-            "area": [0.0, 1000.0] * 2,
+            "area": [0.01, 1000.0] * 2,
             "level": [0.0, 1.0] * 2,
         }
     )

diff --git a/python/ribasim_testmodels/ribasim_testmodels/bucket.py b/python/ribasim_testmodels/ribasim_testmodels/bucket.py
@@ -44,9 +44,9 @@ def bucket_model() -> ribasim.Model:
     # Setup the basins:
     profile = pd.DataFrame(
         data={
-            "node_id": [1, 1, 1],
-            "area": [0.0, 1000.0, 1000.0],
-            "level": [0.0, 0.1, 1.0],
+            "node_id": [1, 1],
+            "area": [1000.0, 1000.0],
+            "level": [0.0, 1.0],
         }
     )
 

diff --git a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py
@@ -56,9 +56,9 @@ def pump_discrete_control_model() -> ribasim.Model:
     # Setup the basins:
     profile = pd.DataFrame(
         data={
-            "node_id": [1, 1, 1, 3, 3, 3],
-            "area": [0.0, 100.0, 100.0] * 2,
-            "level": [0.0, 0.001, 1.0] * 2,
+            "node_id": [1, 1, 3, 3],
+            "area": [100.0, 100.0] * 2,
+            "level": [0.0, 1.0] * 2,
         }
     )
 
@@ -195,9 +195,9 @@ def flow_condition_model():
     # Setup the basins:
     profile = pd.DataFrame(
         data={
-            "node_id": [3, 3, 3],
-            "area": [0.0, 100.0, 100.0],
-            "level": [0.0, 0.001, 1.0],
+            "node_id": [3, 3],
+            "area": [100.0, 100.0],
+            "level": [0.0, 1.0],
         }
     )
 
@@ -331,7 +331,7 @@ def tabulated_rating_curve_control_model() -> ribasim.Model:
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1],
-            "area": [0.0, 1000.0],
+            "area": [0.01, 1000.0],
             "level": [0.0, 1.0],
         }
     )

diff --git a/python/ribasim_testmodels/ribasim_testmodels/equations.py b/python/ribasim_testmodels/ribasim_testmodels/equations.py
@@ -49,7 +49,7 @@ def linear_resistance_model():
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1, 1],
-            "area": [0.0, 100.0, 100.0],
+            "area": [0.01, 100.0, 100.0],
             "level": [0.0, 1.0, 2.0],
         }
     )
@@ -147,7 +147,7 @@ def rating_curve_model():
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1, 1],
-            "area": [0.0, 100.0, 100.0],
+            "area": [0.01, 100.0, 100.0],
             "level": [0.0, 1.0, 2.0],
         }
     )
@@ -258,7 +258,7 @@ def manning_resistance_model():
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1, 1, 3, 3, 3],
-            "area": 2 * [0.0, 100.0, 100.0],
+            "area": 2 * [0.01, 100.0, 100.0],
             "level": 2 * [0.0, 1.0, 2.0],
         }
     )
@@ -366,7 +366,7 @@ def misc_nodes_model():
     profile = pd.DataFrame(
         data={
             "node_id": 3 * [3] + 3 * [5],
-            "area": 2 * [0.0, 100.0, 100.0],
+            "area": 2 * [0.01, 100.0, 100.0],
             "level": 2 * [0.0, 1.0, 2.0],
         }
     )
@@ -500,7 +500,7 @@ def pid_control_equation_model():
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1, 1],
-            "area": [0.0, 100.0, 100.0],
+            "area": [0.01, 100.0, 100.0],
             "level": [0.0, 1.0, 2.0],
         }
     )

diff --git a/python/ribasim_testmodels/ribasim_testmodels/invalid.py b/python/ribasim_testmodels/ribasim_testmodels/invalid.py
@@ -45,7 +45,7 @@ def invalid_qh_model():
     profile = pd.DataFrame(
         data={
             "node_id": [3, 3],
-            "area": [0.0, 1.0],
+            "area": [0.01, 1.0],
             "level": [0.0, 1.0],
         }
     )

diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
@@ -60,6 +60,7 @@ def pid_control_model():
 
     level = np.linspace(0, R, n)
     area = np.pi * level * (2 * R - level)
+    area[0] = 0.01
 
     profile = pd.DataFrame(data={"node_id": n * [2], "level": level, "area": area})
 

diff --git a/python/ribasim_testmodels/ribasim_testmodels/time.py b/python/ribasim_testmodels/ribasim_testmodels/time.py
@@ -50,7 +50,7 @@ def flow_boundary_time_model():
     profile = pd.DataFrame(
         data={
             "node_id": [2, 2],
-            "area": [0.0, 1000.0],
+            "area": [0.01, 1000.0],
             "level": [0.0, 1.0],
         }
     )

diff --git a/python/ribasim_testmodels/ribasim_testmodels/trivial.py b/python/ribasim_testmodels/ribasim_testmodels/trivial.py
@@ -51,7 +51,7 @@ def trivial_model() -> ribasim.Model:
     profile = pd.DataFrame(
         data={
             "node_id": [1, 1],
-            "area": [0.0, 1000.0],
+            "area": [0.01, 1000.0],
             "level": [0.0, 1.0],
         }
     )