Decrease relative tolerance from 1e-5 to 1e-7 (#1874)

With #1819 the state goes further away from 0 with time. That means we
rely more on the relative tolerance than the absolute tolerance, see the
[SciML
docs](https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#What-does-tolerance-mean-and-how-much-error-should-I-expect)
on the topic. Hence I found that reducing the `abstol` had no effect,
but reducing the default `reltol` by two orders of magnitude is enough
to get only negligible leaky Terminals with occasional very brief leaks
of 1e-8 m3/s. It is also enough to minimize the infiltration error to
acceptable levels as can be seen in
#1863 (comment).

~~For HWS this slows down simulation from 9.1s to 12.3s, which isn't too
bad for a two orders of magnitude reduction in relative tolerance.~~
Rechecking HWS performance gives a speedup! 12.4s for 1e-8, 8.1s for
1e-7, 10.2s for 1e-6. So 1e-7 seems optimal. Not sure what went wrong
last time, but this seems pretty consistent. This is with both absolute
and relative set to the same number.

Fixes #1851
Fixes #1863

core/integration_test/hws_integration_test.jl

@@ -17,13 +17,13 @@
 
     @testset "Results values" begin
         @test basin.node_id == basin_bench.node_id
-        @test all(q -> abs(q) < 0.06, basin.level - basin_bench.level)
+        @test all(q -> abs(q) < 0.2, basin.level - basin_bench.level)
     end
 
     timed = @timed Ribasim.run(toml_path)
 
-    # current benchmark is 600s
-    benchmark_runtime = 600
+    # current benchmark in seconds, TeamCity is up to 4x slower than local
+    benchmark_runtime = 32
     performance_diff =
         round((timed.time - benchmark_runtime) / benchmark_runtime * 100; digits = 2)
     if performance_diff < 0.0

core/regression_test/regression_test.jl

@@ -8,7 +8,7 @@
     config = Ribasim.Config(toml_path)
 
     solver_list =
-        ["QNDF", "Rosenbrock23", "TRBDF2", "Rodas5", "KenCarp4", "Tsit5", "ImplicitEuler"]
+        ["QNDF", "Rosenbrock23", "TRBDF2", "Rodas5P", "KenCarp4", "Tsit5", "ImplicitEuler"]
     sparse_on = [true, false]
     autodiff_on = [true, false]
 
@@ -38,13 +38,13 @@
                 # subgrid = Arrow.Table(subgrid_bytes)
 
                 @testset "Results values" begin
-                    @test basin.storage[1] ≈ 1.0
-                    @test basin.level[1] ≈ 0.044711584
-                    @test basin.storage[end] ≈ 16.530443267
+                    @test basin.storage[1] ≈ 1.0f0
+                    @test basin.level[1] ≈ 0.044711584f0
+                    @test basin.storage[end] ≈ 16.530443267f0
                     @test basin.level[end] ≈ 0.181817438
-                    @test flow.flow_rate[1] == basin.outflow_rate[1]
+                    @test flow.flow_rate[1] ≈ basin.outflow_rate[1]
                     @test all(q -> abs(q) < 1e-7, basin.balance_error)
-                    @test all(q -> abs(q) < 0.01, basin.relative_error)
+                    @test all(err -> abs(err) < 0.01, basin.relative_error)
                 end
             end
         end
@@ -105,7 +105,7 @@ end
                 @test all(q -> abs(q) < 1.0, basin.storage - basin_bench.storage)
                 @test all(q -> abs(q) < 0.5, basin.level - basin_bench.level)
                 @test all(q -> abs(q) < 1e-3, basin.balance_error)
-                @test all(q -> abs(q) < 2.5, basin.relative_error)
+                @test all(err -> abs(err) < 2.5, basin.relative_error)
             end
         end
     end
@@ -127,7 +127,7 @@ end
     flow_bench = Arrow.Table(flow_bytes_bench)
     basin_bench = Arrow.Table(basin_bytes_bench)
 
-    # TODO "Rosenbrock23" and "Rodas5" solver are resulting unsolvable gradients
+    # TODO "Rosenbrock23" and "Rodas5P" solver are resulting unsolvable gradients
     solver_list = ["QNDF"]
     sparse_on = [true, false]
     autodiff_on = [true, false]
@@ -167,7 +167,7 @@ end
                 @test basin.node_id == basin_bench.node_id
                 @test all(q -> abs(q) < 100.0, basin.storage - basin_bench.storage)
                 @test all(q -> abs(q) < 0.5, basin.level - basin_bench.level)
-                @test all(q -> abs(q) < 1e-3, basin.balance_error)
+                @test all(err -> abs(err) < 1e-3, basin.balance_error)
             end
         end
     end

core/src/bmi.jl

@@ -48,7 +48,7 @@ function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{F
         p.subgrid.level
     elseif name == "user_demand.demand"
         vec(p.user_demand.demand)
-    elseif name == "user_demand.inflow"
+    elseif name == "user_demand.cumulative_inflow"
         u.user_demand_inflow
     else
         error("Unknown variable $name")

core/src/config.jl

@@ -100,8 +100,8 @@ const nodetypes = collect(keys(nodekinds))
     dtmin::Float64 = 0.0
     dtmax::Union{Float64, Nothing} = nothing
     force_dtmin::Bool = false
-    abstol::Float64 = 1e-6
-    reltol::Float64 = 1e-5
+    abstol::Float64 = 1e-7
+    reltol::Float64 = 1e-7
     water_balance_abstol::Float64 = 1e-3
     water_balance_reltol::Float64 = 1e-2
     maxiters::Int = 1e9

core/test/bmi_test.jl

@@ -11,11 +11,12 @@
     @test BMI.get_end_time(model) ≈ 3.16224e7
     BMI.update(model)
     @test BMI.get_current_time(model) ≈ dt0 atol = 5e-3
-    # cannot go back in time
-    @test_throws ErrorException BMI.update_until(model, dt0 / 2.0)
     @test BMI.get_current_time(model) ≈ dt0 atol = 5e-3
     BMI.update_until(model, 86400.0)
     @test BMI.get_current_time(model) == 86400.0
+    # cannot go back in time
+    @test_throws ErrorException BMI.update_until(model, 3600.0)
+    @test BMI.get_current_time(model) == 86400.0
 end
 
 @testitem "fixed timestepping" begin
@@ -67,7 +68,7 @@ end
         "basin.cumulative_drainage",
         "basin.subgrid_level",
         "user_demand.demand",
-        "user_demand.inflow",
+        "user_demand.cumulative_inflow",
     ]
         value_first = BMI.get_value_ptr(model, name)
         BMI.update_until(model, 86400.0)
@@ -86,7 +87,7 @@ end
     config = Ribasim.Config(toml_path; allocation_use_allocation = false)
     model = Ribasim.Model(config)
     demand = BMI.get_value_ptr(model, "user_demand.demand")
-    inflow = BMI.get_value_ptr(model, "user_demand.inflow")
+    inflow = BMI.get_value_ptr(model, "user_demand.cumulative_inflow")
     # One year in seconds
     year = model.integrator.p.user_demand.demand_itp[2][1].t[2]
     demand_start = 1e-3

core/test/control_test.jl

@@ -247,7 +247,7 @@ end
 
     t_switch = Ribasim.datetime_since(record.time[2], p.starttime)
     flow_table = DataFrame(Ribasim.flow_table(model))
-    @test all(filter(:time => time -> time <= t_switch, flow_table).flow_rate .> 0)
+    @test all(filter(:time => time -> time <= t_switch, flow_table).flow_rate .> -1e-12)
     @test all(
         isapprox.(
             filter(:time => time -> time > t_switch, flow_table).flow_rate,

core/test/docs.toml

@@ -31,8 +31,8 @@ dt = 60.0           # optional, remove for adaptive time stepping
 dtmin = 0.0         # optional, default 0.0
 dtmax = 0.0         # optional, default length of simulation
 force_dtmin = false # optional, default false
-abstol = 1e-6       # optional, default 1e-6
-reltol = 1e-5       # optional, default 1e-5
+abstol = 1e-7       # optional, default 1e-7
+reltol = 1e-7       # optional, default 1e-7
 water_balance_abstol = 1e-3 # optional, default 1e-3
 water_balance_reltol = 1e-2 # optional, default 1e-2
 maxiters = 1e9      # optional, default 1e9

docs/concept/equations.qmd

@@ -157,8 +157,8 @@ For more a more in-depth discussion of numerical computations see [Numerical con
 
 There are many things that can influence the calculations times, for instance:
 
-- [Solver tolerance](https://diffeq.sciml.ai/stable/basics/faq/#What-does-tolerance-mean-and-how-much-error-should-I-expect):
-  By default we use absolute and relative tolerances of `1e-6` and `1e-5` respectively.
+- [Solver tolerance](https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#What-does-tolerance-mean-and-how-much-error-should-I-expect):
+  By default both the absolute and relative tolerance is `1e-7`.
 - [ODE solvers](https://diffeq.sciml.ai/stable/solvers/ode_solve/):
   The `QNDF` method we use is robust to oscillations and massive stiffness, however other solvers should be tried as well.
 - Forcing: Every time new forcing data is injected into the model, it needs to pause.

python/ribasim/ribasim/config.py

@@ -99,9 +99,9 @@ class Solver(ChildModel):
         If a smaller dt than dtmin is needed to meet the set error tolerances, the simulation stops, unless force_dtmin = true
         (Optional, defaults to False)
     abstol : float
-        The absolute tolerance for adaptive timestepping (Optional, defaults to 1e-6)
+        The absolute tolerance for adaptive timestepping (Optional, defaults to 1e-7)
     reltol : float
-        The relative tolerance for adaptive timestepping (Optional, defaults to 1e-5)
+        The relative tolerance for adaptive timestepping (Optional, defaults to 1e-7)
     maxiters : int
         The total number of linear iterations over the whole simulation. (Defaults to 1e9, only needs to be increased for extremely long simulations)
     sparse : bool

python/ribasim_api/tests/test_bmi.py

@@ -168,7 +168,7 @@ def test_get_value_ptr_allocation(libribasim, user_demand, tmp_path):
     libribasim.update_until(60.0)
 
     # Realized inflow
-    actual_inflow = libribasim.get_value_ptr("user_demand.inflow")
+    actual_inflow = libribasim.get_value_ptr("user_demand.cumulative_inflow")
     expected_inflow = np.array([1e-4 * 60.0, 0.0, 0.0])
     assert_array_almost_equal(actual_inflow, expected_inflow)