Remove adaptive solver config option (#1170)

Instead, the presence of `dt` will determine if it is adaptive or not.

The only feature we lose with this is that we cannot set the initial
timestep for an adaptive run, which to my knowledge wasn't used.

# Upgrade information

In the TOML, remove the `adaptive` key in the `[solver]` section.
Provide a `dt` key only if you want a fixed timestep.

---------

Co-authored-by: Jingru Feng <[email protected]>

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, convert_saveat
+export algorithm, snake_case, input_path, results_path, convert_saveat, convert_dt
 
 const schemas =
     getfield.(
@@ -78,7 +78,6 @@ const nodetypes = collect(keys(nodekinds))
 @option struct Solver <: TableOption
     algorithm::String = "QNDF"
     saveat::Float64 = 86400.0
-    adaptive::Bool = true
     dt::Union{Float64, Nothing} = nothing
     dtmin::Float64 = 0.0
     dtmax::Union{Float64, Nothing} = nothing
@@ -248,4 +247,23 @@ function convert_saveat(saveat::Float64, t_end::Float64)::Union{Float64, Vector{
     end
 end
 
+"Convert the dt from our Config to SciML stepsize control arguments"
+function convert_dt(dt::Union{Float64, Nothing})::Tuple{Bool, Float64}
+    # In SciML dt represents the initial timestep if adaptive is true.
+    # We don't support setting the initial timestep, so we don't need the adaptive flag.
+    # The solver will give a clear error message if the algorithm is not adaptive.
+    if isnothing(dt)
+        # adaptive step size
+        adaptive = true
+        dt = 0.0
+    elseif 0 < dt < Inf
+        # fixed step size
+        adaptive = false
+    else
+        @error "Invalid dt" dt
+        error("Invalid dt")
+    end
+    adaptive, dt
+end
+
 end  # module

core/src/model.jl

@@ -106,6 +106,7 @@ function Model(config::Config)::Model
     t0 = zero(t_end)
     timespan = (t0, t_end)
     saveat = convert_saveat(config.solver.saveat, t_end)
+    adaptive, dt = convert_dt(config.solver.dt)
 
     jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing
     RHS = ODEFunction(water_balance!; jac_prototype)
@@ -132,8 +133,8 @@ function Model(config::Config)::Model
         tstops,
         isoutofdomain = (u, p, t) -> any(<(0), u.storage),
         saveat,
-        config.solver.adaptive,
-        dt = something(config.solver.dt, t0),
+        adaptive,
+        dt,
         config.solver.dtmin,
         dtmax = something(config.solver.dtmax, t_end),
         config.solver.force_dtmin,

core/test/bmi_test.jl

@@ -23,14 +23,9 @@ end
 
     toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml")
     dt = 10.0
-    config = Ribasim.Config(
-        toml_path;
-        solver_algorithm = "ImplicitEuler",
-        solver_adaptive = false,
-        solver_dt = dt,
-    )
+    config = Ribasim.Config(toml_path; solver_algorithm = "ImplicitEuler", solver_dt = dt)
     @test config.solver.algorithm == "ImplicitEuler"
-    @test !config.solver.adaptive
+    @test config.solver.dt === dt
     model = Ribasim.Model(config)
 
     @test BMI.get_time_step(model) == dt

core/test/config_test.jl

@@ -32,21 +32,20 @@
     @testset "docs" begin
         config = Ribasim.Config(normpath(@__DIR__, "docs.toml"))
         @test config isa Ribasim.Config
-        @test config.solver.adaptive
+        @test config.solver.autodiff
     end
 end
 
 @testitem "Solver" begin
     using OrdinaryDiffEq: alg_autodiff, AutoFiniteDiff, AutoForwardDiff
-    using Ribasim: convert_saveat, Solver, algorithm
+    using Ribasim: convert_saveat, convert_dt, Solver, algorithm
 
     solver = Solver()
     @test solver.algorithm == "QNDF"
     Solver(;
         algorithm = "Rosenbrock23",
         autodiff = true,
         saveat = 3600.0,
-        adaptive = true,
         dt = 0,
         abstol = 1e-5,
         reltol = 1e-4,
@@ -72,6 +71,10 @@ 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)
+
+    @test convert_dt(nothing) == (true, 0.0)
+    @test convert_dt(360.0) == (false, 360.0)
+    @test_throws ErrorException convert_dt(0.0)
 end
 
 @testitem "snake_case" begin

core/test/docs.toml

@@ -21,8 +21,7 @@ objective_type = "linear_absolute" # optional, default "linear_absolute"
 [solver]
 algorithm = "QNDF"  # optional, default "QNDF"
 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
+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

core/test/equations_test.jl

@@ -168,10 +168,12 @@ end
 @testitem "MiscellaneousNodes" begin
     using PreallocationTools: get_tmp
     using SciMLBase: successful_retcode
+    using Ribasim: timesteps, get_storages_and_levels
 
     toml_path = normpath(@__DIR__, "../../generated_testmodels/misc_nodes/ribasim.toml")
     @test ispath(toml_path)
-    model = Ribasim.run(toml_path)
+    config = Ribasim.Config(toml_path)
+    model = Ribasim.run(config)
     @test successful_retcode(model)
     p = model.integrator.p
     (; flow_boundary, fractional_flow, pump) = p
@@ -181,9 +183,12 @@ end
     q_pump = pump_flow_rate[1]
     frac = fractional_flow.fraction[1]
 
-    storage_both = Ribasim.get_storages_and_levels(model).storage
-    t = Ribasim.timesteps(model)
+    storage_both = get_storages_and_levels(model).storage
+    t = timesteps(model)
+    tspan = model.integrator.sol.prob.tspan
 
+    @test config.solver.dt === model.integrator.dt
+    @test t ≈ range(tspan...; step = config.solver.saveat)
     @test storage_both[1, :] ≈ @. storage_both[1, 1] + t * (frac * q_boundary - q_pump)
     @test storage_both[2, :] ≈ @. storage_both[2, 1] + t * q_pump
 end

docs/core/usage.qmd

@@ -31,19 +31,17 @@ Ribasim expects the GeoPackage database [database.gpkg](usage.qmd#sec-geopackage
 ## Solver settings {#sec-solver-settings}
 
 The solver section in the configuration file is entirely optional, since we aim to use defaults that will generally work well.
-Common reasons to modify the solver settings are to adjust the calculation or result stepsizes: `adaptive`, `dt`, and `saveat`.
+Common reasons to modify the solver settings are to adjust the calculation or result stepsizes: `dt`, and `saveat`.
 If your model does not converge, or your performance is lower than expected, it can help to adjust other solver settings as well.
 
 The default solver `algorithm = "QNDF"`, which is a multistep method similar to Matlab's `ode15s` [@shampine1997matlab].
-It is an implicit method that supports the default `adaptive = true` timestepping.
+It is an implicit method that supports the default adaptive timestepping.
 The full list of available solvers is: `QNDF`, `Rosenbrock23`, `TRBDF2`, `Rodas5`, `KenCarp4`, `Tsit5`, `RK4`, `ImplicitEuler`, `Euler`.
 Information on the solver algorithms can be found on the [ODE solvers page](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).
 
-The `dt` controls the stepsize.
-When `adaptive = true`, `dt` only applies to the initial stepsize, and by default it is automatically determined.
-When `adaptive = false` a suitable `dt` must always be provided.
-The value is in seconds, so `dt = 3600.0` corresponds to an hourly timestep.
-When `adaptive = true`, `dtmin` and `dtmax` control the minimum and maximum allowed `dt`.
+By default Ribasim uses adaptive timestepping, though not all algorithms support adaptive timestepping.
+To use fixed timesteps, provide a timestep size in seconds; `dt = 3600.0` corresponds to an hourly timestep.
+With adaptive timestepping, `dtmin` and `dtmax` control the minimum and maximum allowed `dt`.
 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.
 

python/ribasim/ribasim/config.py

@@ -50,7 +50,6 @@ class Results(ChildModel):
 class Solver(ChildModel):
     algorithm: str = "QNDF"
     saveat: float = 86400.0
-    adaptive: bool = True
     dt: float | None = None
     dtmin: float | None = None
     dtmax: float | None = None

python/ribasim_testmodels/ribasim_testmodels/equations.py

@@ -384,7 +384,6 @@ def misc_nodes_model():
 
     # Setup solver:
     solver = ribasim.Solver(
-        adaptive=False,
         dt=24 * 24 * 60,
         algorithm="Euler",
     )