lower reltol from 1e-3 to 1e-5 (#773)

@harm-nomden-sweco shared a case where a pump is emptying a basin with
an inflow smaller than the pump capacity. We know this case leads to
small timesteps, but it appeared this also led to large errors in the
water balance. This lowering of the reltol is enough to fix that issue.
Since such a case can be common, it is good to have a safe default.

This also adds `using Test` to all test files. And `test/time.jl` needed
some refactoring since there were a few sub-millisecond timesteps,
whereas we round our time to milliseconds.

core/src/config.jl

@@ -84,7 +84,7 @@ const nodetypes = collect(keys(nodekinds))
     dtmax::Union{Float64, Nothing} = nothing
     force_dtmin::Bool = false
     abstol::Float64 = 1e-6
-    reltol::Float64 = 1e-3
+    reltol::Float64 = 1e-5
     maxiters::Int = 1e9
     sparse::Bool = true
     autodiff::Bool = true

core/test/allocation.jl

@@ -1,3 +1,4 @@
+using Test
 import Ribasim
 import JuMP
 using SQLite

core/test/bmi.jl

@@ -7,14 +7,15 @@ toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml")
 @testset "adaptive timestepping" begin
     model = BMI.initialize(Ribasim.Model, toml_path)
     @test BMI.get_time_units(model) == "s"
-    dt0 = 0.011371289f0
+    dt0 = 0.0001269439f0
     @test BMI.get_time_step(model) ≈ dt0 atol = 5e-3
     @test BMI.get_start_time(model) === 0.0
     @test BMI.get_current_time(model) === 0.0
     @test BMI.get_end_time(model) ≈ 3.16224e7
     BMI.update(model)
     @test BMI.get_current_time(model) ≈ dt0 atol = 5e-3
-    @test_throws ErrorException BMI.update_until(model, 0.005)
+    # 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

core/test/control.jl

@@ -1,3 +1,4 @@
+using Test
 import Ribasim
 using Dates: Date
 

core/test/docs.jl

@@ -1,4 +1,7 @@
-using Test, Documenter, Ribasim
+using Test
+using Documenter
+using Ribasim
+
 DocMeta.setdocmeta!(Ribasim, :DocTestSetup, :(using Ribasim); recursive = true)
 
 @testset "Doctests" begin

core/test/docs.toml

@@ -25,7 +25,7 @@ dtmin = 0.0         # optional, default automatically determined
 dtmax = 0.0         # optional, default length of simulation
 force_dtmin = false # optional, default false
 abstol = 1e-6       # optional, default 1e-6
-reltol = 1e-3       # optional, default 1e-3
+reltol = 1e-5       # optional, default 1e-5
 maxiters = 1e9      # optional, default 1e9
 sparse = true       # optional, default true
 autodiff = true     # optional, default true

core/test/run_models.jl

@@ -44,7 +44,7 @@ end
     @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip =
         Sys.isapple() atol = 1.5
 
-    @test length(logger.logs) == 7
+    @test length(logger.logs) == 8
     @test logger.logs[1].level == Debug
     @test logger.logs[1].message == "Read database into memory."
 
@@ -74,7 +74,7 @@ end
     @test model isa Ribasim.Model
     @test successful_retcode(model)
     @test length(model.integrator.p.basin.precipitation) == 4
-    @test model.integrator.sol.u[end] ≈ Float32[469.8923, 469.89038, 410.71472, 1427.4194] skip =
+    @test model.integrator.sol.u[end] ≈ Float32[472.02444, 472.02252, 367.6387, 1427.981] skip =
         Sys.isapple()
 end
 
@@ -198,8 +198,8 @@ end
     @test only(model.integrator.sol(0day)) == 1000.0
     # constant user withdraws to 0.9m/900m3
     @test only(model.integrator.sol(150day)) ≈ 900 atol = 5
-    # dynamic user withdraws to 0.5m/500m3
-    @test only(model.integrator.sol(180day)) ≈ 500 atol = 1
+    # dynamic user withdraws to 0.5m/509m3
+    @test only(model.integrator.sol(180day)) ≈ 509 atol = 1
 end
 
 @testset "ManningResistance" begin

core/test/time.jl

@@ -1,7 +1,7 @@
 using Ribasim
-using SQLite
 using Dates
 using DataFrames: DataFrame
+using Test
 
 @testset "Time dependent flow boundary" begin
     toml_path =
@@ -11,19 +11,22 @@ using DataFrames: DataFrame
 
     flow = DataFrame(Ribasim.flow_table(model))
     # only from March to September the FlowBoundary varies
-    sin_timestamps = 3 .<= month.(unique(flow.time)) .< 10
+    is_summer(t::DateTime) = 3 <= month(t) < 10
 
-    flow_added_1 = filter(
-        [:from_node_id, :to_node_id] => (from, to) -> from === 1 && to === 1,
-        flow,
-    ).flow[sin_timestamps]
+    flow_added_1 =
+        filter(
+            [:time, :from_node_id, :to_node_id] =>
+                (t, from, to) -> is_summer(t) && from === 1 && to === 1,
+            flow,
+        ).flow
     flow_1_to_2 = filter(
-        [:from_node_id, :to_node_id] => (from, to) -> from === 1 && to === 2,
+        [:time, :from_node_id, :to_node_id] =>
+            (t, from, to) -> is_summer(t) && from === 1 && to === 2,
         flow,
-    ).flow[sin_timestamps]
-    @test flow_added_1 == flow_1_to_2
+    )
+    @test flow_added_1 == flow_1_to_2.flow
 
-    t = model.saved_flow.t[sin_timestamps]
+    t = Ribasim.seconds_since.(flow_1_to_2.time, model.config.starttime)
     flow_expected = @. 1 + sin(0.5 * π * (t - t[1]) / (t[end] - t[1]))^2
     # some difference is expected since the modeled flow is for the period up to t
     @test isapprox(flow_added_1, flow_expected, rtol = 0.005)

core/test/validation.jl

@@ -1,3 +1,4 @@
+using Test
 using Ribasim
 using Graphs: DiGraph, add_edge!
 using Dictionaries: Indices

docs/schema/Config.schema.json

@@ -51,7 +51,7 @@
                 "dtmax": null,
                 "force_dtmin": false,
                 "abstol": 1.0e-6,
-                "reltol": 0.001,
+                "reltol": 1.0e-5,
                 "maxiters": 1000000000,
                 "sparse": true,
                 "autodiff": true

docs/schema/Solver.schema.json

@@ -80,7 +80,7 @@
         "reltol": {
             "format": "double",
             "type": "number",
-            "default": 0.001
+            "default": 1.0e-5
         },
         "maxiters": {
             "format": "default",

python/ribasim/ribasim/config.py

@@ -23,7 +23,7 @@ class Solver(BaseModel):
     dtmax: Optional[float] = None
     force_dtmin: bool = False
     abstol: float = 1e-06
-    reltol: float = 0.001
+    reltol: float = 1e-05
     maxiters: int = 1000000000
     sparse: bool = True
     autodiff: bool = True
@@ -128,7 +128,7 @@ class Config(BaseModel):
                 "dtmax": None,
                 "force_dtmin": False,
                 "abstol": 1e-06,
-                "reltol": 0.001,
+                "reltol": 1e-05,
                 "maxiters": 1000000000,
                 "sparse": True,
                 "autodiff": True,