Show convergence bottleneck nodes to users

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.0"
 manifest_format = "2.0"
-project_hash = "045162196d64ccb59bdab2af0422fde0b0e83c77"
+project_hash = "84ca8e19d9af141a7622a8feb4f64903c913b47b"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "016833eb52ba2d6bea9fcb50ca295980e728ee24"
@@ -1525,6 +1525,12 @@ deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
 version = "7.2.1+1"
 
+[[deps.Suppressor]]
+deps = ["Logging"]
+git-tree-sha1 = "9143c41bd539a8885c79728b9dedb0ce47dc9819"
+uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+version = "0.2.7"
+
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "MacroTools", "RuntimeGeneratedFunctions", "StaticArraysCore"]
 git-tree-sha1 = "4b7f4c80449d8baae8857d55535033981862619c"

Project.toml

@@ -49,6 +49,7 @@ SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"

core/Project.toml

@@ -93,6 +93,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IOCapture = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

core/src/main.jl

@@ -46,6 +46,30 @@ function main(toml_path::AbstractString)::Cint
                 t = datetime_since(model.integrator.t, starttime)
                 retcode = model.integrator.sol.retcode
                 @error "The model exited at model time $t with return code $retcode.\nSee https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#retcodes"
+
+                (; cache, p, u) = model.integrator
+                (; basin) = p
+                # Indicate convergence bottlenecks if possible with the current algorithm
+                if hasproperty(cache, :atmp)
+                    storage_error = @. abs(cache.atmp.storage / u.storage)
+                    perm = sortperm(storage_error; rev = true)
+                    println(
+                        "The following basins were identified as convergence bottlenecks (in descending order of severity):",
+                    )
+                    for i in perm
+                        node_id = basin.node_id.values[i]
+                        error = storage_error[i]
+                        if error < config.solver.reltol
+                            break
+                        end
+                        println("$node_id, error = $error")
+                    end
+                else
+                    algorithm = model.config.solver_algorithm
+                    println(
+                        "The current algorithm ($algorithm) does not support indicating convergence bottlenecks, for that try a different one.",
+                    )
+                end
                 return 1
             end
         end

core/test/validation_test.jl

@@ -458,3 +458,19 @@ end
         dt,
     )
 end
+
+@testitem "Convergence bottleneck" begin
+    using Suppressor: @capture_out
+    toml_path =
+        normpath(@__DIR__, "../../generated_testmodels/invalid_unstable/ribasim.toml")
+    @test ispath(toml_path)
+    output = @capture_out Ribasim.main(toml_path)
+    output = split(output, "\n")
+    @test startswith(
+        output[1],
+        "The following basins were identified as convergence bottlenecks",
+    )
+    @test startswith(output[2], "Basin #1")
+    @test startswith(output[3], "Basin #21")
+    @test startswith(output[4], "Basin #41")
+end

python/ribasim_testmodels/ribasim_testmodels/invalid.py

@@ -222,22 +222,26 @@ def invalid_edge_types_model() -> Model:
 
 
 def invalid_unstable_model() -> Model:
-    """Model with an extremely quickly emptying basin."""
+    """Model with several extremely quickly emptying basins."""
 
     model = Model(
         starttime="2020-01-01",
         endtime="2021-01-01",
         crs="EPSG:28992",
-        solver=Solver(dtmin=1.0),
-    )
-
-    model.basin.add(
-        Node(1, Point(0, 0)),
-        [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[1.0])],
-    )
-    model.pump.add(Node(2, Point(0, 1)), [pump.Static(flow_rate=[1e15])])
-    model.terminal.add(Node(3, Point(0, 2)))
-
-    model.edge.add(model.basin[1], model.pump[2])
-    model.edge.add(model.pump[2], model.terminal[3])
+        solver=Solver(dtmin=60.0),
+    )
+    id_shift = 10
+    for i in range(6):
+        model.basin.add(
+            Node(1 + id_shift * i, Point(i, 0)),
+            [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[1.0])],
+        )
+        flow_rate = 1.0 if (i % 2 == 0) else 1e15
+        model.pump.add(
+            Node(2 + id_shift * i, Point(i, 1)), [pump.Static(flow_rate=[flow_rate])]
+        )
+        model.terminal.add(Node(3 + id_shift * i, Point(i, 2)))
+
+        model.edge.add(model.basin[1 + id_shift * i], model.pump[2 + id_shift * i])
+        model.edge.add(model.pump[2 + id_shift * i], model.terminal[3 + id_shift * i])
     return model