Merge branch 'main' into mean_flow_output

core/src/config.jl

@@ -16,7 +16,7 @@ using ..Ribasim: Ribasim, isnode, nodetype
 using OrdinaryDiffEq
 
 export Config, Solver, Results, Logging, Toml
-export algorithm, snake_case, input_path, results_path
+export algorithm, snake_case, input_path, results_path, convert_saveat
 
 const schemas =
     getfield.(
@@ -77,7 +77,7 @@ const nodetypes = collect(keys(nodekinds))
 
 @option struct Solver <: TableOption
     algorithm::String = "QNDF"
-    saveat::Union{Float64, Vector{Float64}} = Float64[]
+    saveat::Float64 = 86400.0
     adaptive::Bool = true
     dt::Union{Float64, Nothing} = nothing
     dtmin::Float64 = 0.0
@@ -166,11 +166,6 @@ function Configurations.from_dict(::Type{Logging}, ::Type{LogLevel}, level::Abst
     )
 end
 
-# [] in TOML is parsed as a Vector{Union{}}
-function Configurations.from_dict(::Type{Solver}, t::Type, saveat::Vector{Union{}})
-    return Float64[]
-end
-
 # TODO Use with proper alignment
 function Base.show(io::IO, c::Config)
     println(io, "Ribasim Config")
@@ -236,4 +231,21 @@ function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm
     end
 end
 
+"Convert the saveat Float64 from our Config to SciML's saveat"
+function convert_saveat(saveat::Float64, t_end::Float64)::Union{Float64, Vector{Float64}}
+    if iszero(saveat)
+        # every step
+        Float64[]
+    elseif saveat == Inf
+        # only the start and end
+        [0.0, t_end]
+    elseif isfinite(saveat)
+        # every saveat seconds
+        saveat
+    else
+        @error "Invalid saveat" saveat
+        error("Invalid saveat")
+    end
+end
+
 end  # module

core/src/model.jl

@@ -105,6 +105,7 @@ function Model(config::Config)::Model
     @assert eps(t_end) < 3600 "Simulation time too long"
     t0 = zero(t_end)
     timespan = (t0, t_end)
+    saveat = convert_saveat(config.solver.saveat, t_end)
 
     jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing
     RHS = ODEFunction(water_balance!; jac_prototype)
@@ -114,7 +115,7 @@ function Model(config::Config)::Model
     end
     @debug "Setup ODEProblem."
 
-    callback, saved = create_callbacks(parameters, config; config.solver.saveat)
+    callback, saved = create_callbacks(parameters, config; saveat)
     @debug "Created callbacks."
 
     # Initialize the integrator, providing all solver options as described in
@@ -130,7 +131,7 @@ function Model(config::Config)::Model
         callback,
         tstops,
         isoutofdomain = (u, p, t) -> any(<(0), u.storage),
-        config.solver.saveat,
+        saveat,
         config.solver.adaptive,
         dt = something(config.solver.dt, t0),
         config.solver.dtmin,

core/test/config_test.jl

@@ -7,7 +7,7 @@
         config = Ribasim.Config(normpath(@__DIR__, "data", "config_test.toml"))
         @test config isa Ribasim.Config
         @test config.endtime > config.starttime
-        @test config.solver == Ribasim.Solver(; saveat = 86400.0)
+        @test config.solver == Ribasim.Solver(; saveat = 3600.0)
         @test config.results.compression
         @test config.results.compression_level == 6
     end
@@ -38,10 +38,11 @@ end
 
 @testitem "Solver" begin
     using OrdinaryDiffEq: alg_autodiff, AutoFiniteDiff, AutoForwardDiff
+    using Ribasim: convert_saveat, Solver, algorithm
 
-    solver = Ribasim.Solver()
+    solver = Solver()
     @test solver.algorithm == "QNDF"
-    Ribasim.Solver(;
+    Solver(;
         algorithm = "Rosenbrock23",
         autodiff = true,
         saveat = 3600.0,
@@ -51,21 +52,26 @@ end
         reltol = 1e-4,
         maxiters = 1e5,
     )
-    Ribasim.Solver(; algorithm = "DoesntExist")
-    @test_throws InexactError Ribasim.Solver(autodiff = 2)
-    @test_throws "algorithm DoesntExist not supported" Ribasim.algorithm(
-        Ribasim.Solver(; algorithm = "DoesntExist"),
+    Solver(; algorithm = "DoesntExist")
+    @test_throws InexactError Solver(autodiff = 2)
+    @test_throws "algorithm DoesntExist not supported" algorithm(
+        Solver(; algorithm = "DoesntExist"),
     )
-    @test alg_autodiff(
-        Ribasim.algorithm(Ribasim.Solver(; algorithm = "QNDF", autodiff = true)),
-    ) == AutoForwardDiff()
-    @test alg_autodiff(
-        Ribasim.algorithm(Ribasim.Solver(; algorithm = "QNDF", autodiff = false)),
-    ) == AutoFiniteDiff()
-    @test alg_autodiff(Ribasim.algorithm(Ribasim.Solver(; algorithm = "QNDF"))) ==
+    @test alg_autodiff(algorithm(Solver(; algorithm = "QNDF", autodiff = true))) ==
           AutoForwardDiff()
+    @test alg_autodiff(algorithm(Solver(; algorithm = "QNDF", autodiff = false))) ==
+          AutoFiniteDiff()
+    @test alg_autodiff(algorithm(Solver(; algorithm = "QNDF"))) == AutoForwardDiff()
     # autodiff is not a kwargs for explicit algorithms, but we use try-catch to bypass
-    Ribasim.algorithm(Ribasim.Solver(; algorithm = "Euler", autodiff = true))
+    algorithm(Solver(; algorithm = "Euler", autodiff = true))
+
+    t_end = 100.0
+    @test convert_saveat(0.0, t_end) == Float64[]
+    @test convert_saveat(60.0, t_end) == 60.0
+    @test convert_saveat(Inf, t_end) == [0.0, t_end]
+    @test convert_saveat(Inf, t_end) == [0.0, t_end]
+    @test_throws ErrorException convert_saveat(-Inf, t_end)
+    @test_throws ErrorException convert_saveat(NaN, t_end)
 end
 
 @testitem "snake_case" begin

core/test/data/config_test.toml

@@ -11,4 +11,4 @@ database = "database.gpkg"
 time = "basin/time.arrow"
 
 [solver]
-saveat = 86400
+saveat = 3600

core/test/docs.toml

@@ -20,7 +20,7 @@ objective_type = "linear_absolute" # optional, default "linear_absolute"
 
 [solver]
 algorithm = "QNDF"  # optional, default "QNDF"
-saveat = []         # optional, default [], which saves every timestep
+saveat = 86400      # optional, default 86400, 0 saves every timestep, inf saves only at start- and endtime
 adaptive = true     # optional, default true
 dt = 0.0            # optional when adaptive = true, default automatically determined
 dtmin = 0.0         # optional, default 0.0

docs/core/usage.qmd

@@ -47,9 +47,10 @@ When `adaptive = true`, `dtmin` and `dtmax` control the minimum and maximum allo
 If a smaller `dt` than `dtmin` is needed to meet the set error tolerances, the simulation stops, unless `force_dtmin` is set to `true`.
 `force_dtmin` is off by default to ensure an accurate solution.
 
-By default the calculation and result stepsize are the same, with `saveat = []`, which will save every timestep.
-`saveat` can be a number, which is the saving interval in seconds, or it can be a list of numbers, which are the times in seconds since start that are saved.
-For instance, `saveat = 86400.0` will save results after every day that passed.
+The default result stepsize, `saveat = 86400` will save results after every day that passed.
+The calculation and result stepsize need not be the same.
+If you wish to save every calculation step, set `saveat = 0`.
+If you wish to not save any intermediate steps, set `saveat = inf`.
 
 The Jacobian matrix provides information about the local sensitivity of the model with respect to changes in the states.
 For implicit solvers it must be calculated often, which can be expensive to do.

python/ribasim/ribasim/config.py

@@ -49,7 +49,7 @@ class Results(ChildModel):
 
 class Solver(ChildModel):
     algorithm: str = "QNDF"
-    saveat: float | list[float] = []
+    saveat: float = 86400.0
     adaptive: bool = True
     dt: float | None = None
     dtmin: float | None = None

python/ribasim/tests/test_model.py

@@ -17,13 +17,16 @@ def test_repr(basic):
 def test_solver():
     solver = Solver()
     assert solver.algorithm == "QNDF"  # default
-    assert solver.saveat == []
+    assert solver.saveat == 86400.0
 
     solver = Solver(saveat=3600.0)
     assert solver.saveat == 3600.0
 
-    solver = Solver(saveat=[3600.0, 7200.0])
-    assert solver.saveat == [3600.0, 7200.0]
+    solver = Solver(saveat=float("inf"))
+    assert solver.saveat == float("inf")
+
+    solver = Solver(saveat=0)
+    assert solver.saveat == 0
 
     with pytest.raises(ValidationError):
         Solver(saveat="a")

python/ribasim_testmodels/ribasim_testmodels/basic.py

@@ -481,6 +481,7 @@ def outlet_model():
         level_boundary=level_boundary,
         starttime="2020-01-01 00:00:00",
         endtime="2021-01-01 00:00:00",
+        solver=ribasim.Solver(saveat=0),
     )
 
     return model