Merge branch 'main' into allocation

core/src/bmi.jl

@@ -105,7 +105,8 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
     t_end = seconds_since(config.endtime, config.starttime)
     # for Float32 this method allows max ~1000 year simulations without accuracy issues
     @assert eps(t_end) < 3600 "Simulation time too long"
-    timespan = (zero(t_end), t_end)
+    t0 = zero(t_end)
+    timespan = (t0, t_end)
 
     jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing
     RHS = ODEFunction(water_balance!; jac_prototype)
@@ -118,6 +119,10 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
     callback, saved_flow = create_callbacks(parameters; config.solver.saveat)
     @debug "Created callbacks."
 
+    # Initialize the integrator, providing all solver options as described in
+    # https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/
+    # Not all keyword arguments (e.g. `dt`) support `nothing`, in which case we follow
+    # https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.57.0/src/solve.jl#L10
     @timeit_debug to "Setup integrator" integrator = init(
         prob,
         alg;
@@ -126,9 +131,13 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
         progress_steps = 100,
         callback,
         tstops,
+        isoutofdomain = (u, p, t) -> any(<(0), u.storage),
         config.solver.saveat,
         config.solver.adaptive,
-        config.solver.dt,
+        dt = something(config.solver.dt, t0),
+        config.solver.dtmin,
+        dtmax = something(config.solver.dtmax, t_end),
+        config.solver.force_dtmin,
         config.solver.abstol,
         config.solver.reltol,
         config.solver.maxiters,

core/src/config.jl

@@ -8,7 +8,7 @@ Ribasim.config is a submodule mainly to avoid name clashes between the configura
 """
 module config
 
-using Configurations: Configurations, Maybe, @option, from_toml, @type_alias
+using Configurations: Configurations, @option, from_toml, @type_alias
 using DataStructures: DefaultDict
 using Dates
 using Logging: LogLevel, Debug, Info, Warn, Error
@@ -41,11 +41,14 @@ end
 snake_case(sym::Symbol)::Symbol = Symbol(snake_case(String(sym)))
 
 """
-Add fieldnames with Maybe{String} type to struct expression. Requires @option use before it.
+Add fieldnames with Union{String, Nothing} type to struct expression. Requires @option use before it.
 """
 macro addfields(typ::Expr, fieldnames)
     for fieldname in fieldnames
-        push!(typ.args[3].args, Expr(:(=), Expr(:(::), fieldname, Maybe{String}), nothing))
+        push!(
+            typ.args[3].args,
+            Expr(:(=), Expr(:(::), fieldname, Union{String, Nothing}), nothing),
+        )
     end
     return esc(typ)
 end
@@ -76,7 +79,10 @@ const nodetypes = collect(keys(nodekinds))
     algorithm::String = "QNDF"
     saveat::Union{Float64, Vector{Float64}} = Float64[]
     adaptive::Bool = true
-    dt::Float64 = 0.0
+    dt::Union{Float64, Nothing} = nothing
+    dtmin::Union{Float64, Nothing} = nothing
+    dtmax::Union{Float64, Nothing} = nothing
+    force_dtmin::Bool = false
     abstol::Float64 = 1e-6
     reltol::Float64 = 1e-3
     maxiters::Int = 1e9
@@ -106,7 +112,7 @@ end
     basin::String = "output/basin.arrow"
     flow::String = "output/flow.arrow"
     control::String = "output/control.arrow"
-    outstate::Maybe{String}
+    outstate::Union{String, Nothing} = nothing
     compression::Compression = "zstd"
     compression_level::Int = 6
 end

core/src/solve.jl

@@ -1169,14 +1169,6 @@ function water_balance!(
     # PID control (changes the du of PID controlled basins)
     continuous_control!(u, du, pid_control, p, integral, t)
 
-    # Negative storage must not decrease, based on Shampine's et. al. advice
-    # https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain
-    for i in eachindex(u.storage)
-        if u.storage[i] < 0
-            du.storage[i] = max(du.storage[i], 0.0)
-        end
-    end
-
     return nothing
 end
 

core/test/docs.toml

@@ -24,15 +24,18 @@ compression_level = 6            # optional, default 6
 time = "basin/time.arrow"
 
 [solver]
-algorithm = "QNDF" # optional, default "QNDF"
-saveat = []        # optional, default [], which saves every timestep
-adaptive = true    # optional, default true
-dt = 0.0           # optional, default 0.0
-abstol = 1e-6      # optional, default 1e-6
-reltol = 1e-3      # optional, default 1e-3
-maxiters = 1e9     # optional, default 1e9
-sparse = true      # optional, default true
-autodiff = true    # optional, default true
+algorithm = "QNDF"  # optional, default "QNDF"
+saveat = []         # optional, default [], which saves every timestep
+adaptive = true     # optional, default true
+dt = 0.0            # optional when adaptive = true, default automatically determined
+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
+maxiters = 1e9      # optional, default 1e9
+sparse = true       # optional, default true
+autodiff = true     # optional, default true
 
 [logging]
 # defines the logging level of Ribasim

core/test/run_models.jl

@@ -38,6 +38,7 @@ end
     @test all(isconcretetype, fieldtypes(typeof(p)))
 
     @test successful_retcode(model)
+    @test allunique(model.integrator.sol.t)
     @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip =
         Sys.isapple() atol = 1.5
 

docs/core/usage.qmd

@@ -40,9 +40,13 @@ The full list of available solvers is: `QNDF`, `Rosenbrock23`, `TRBDF2`, `Rodas5
 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` this only applies to the initial stepsize, and the default of 0.0 means that it is automatically determined.
+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 these depend on the problem and algorithm.
+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 output 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.

docs/schema/Config.schema.json

@@ -64,7 +64,10 @@
                 "saveat": [
                 ],
                 "adaptive": true,
-                "dt": 0,
+                "dt": null,
+                "dtmin": null,
+                "dtmax": null,
+                "force_dtmin": false,
                 "abstol": 1.0e-6,
                 "reltol": 0.001,
                 "maxiters": 1000000000,

docs/schema/Solver.schema.json

@@ -1,10 +1,32 @@
 {
     "$schema": "https://json-schema.org/draft/2020-12/schema",
     "properties": {
-        "reltol": {
-            "format": "double",
-            "default": 0.001,
-            "type": "number"
+        "autodiff": {
+            "format": "default",
+            "default": true,
+            "type": "boolean"
+        },
+        "adaptive": {
+            "format": "default",
+            "default": true,
+            "type": "boolean"
+        },
+        "force_dtmin": {
+            "format": "default",
+            "default": false,
+            "type": "boolean"
+        },
+        "dt": {
+            "format": "default",
+            "anyOf": [
+                {
+                    "type": "null"
+                },
+                {
+                    "type": "number"
+                }
+            ],
+            "default": null
         },
         "saveat": {
             "format": "default",
@@ -27,12 +49,19 @@
             "default": 1000000000,
             "type": "integer"
         },
-        "autodiff": {
+        "dtmin": {
             "format": "default",
-            "default": true,
-            "type": "boolean"
+            "anyOf": [
+                {
+                    "type": "null"
+                },
+                {
+                    "type": "number"
+                }
+            ],
+            "default": null
         },
-        "adaptive": {
+        "sparse": {
             "format": "default",
             "default": true,
             "type": "boolean"
@@ -47,22 +76,29 @@
             "default": 1.0e-6,
             "type": "number"
         },
-        "dt": {
+        "reltol": {
             "format": "double",
-            "default": 0,
+            "default": 0.001,
             "type": "number"
         },
-        "sparse": {
+        "dtmax": {
             "format": "default",
-            "default": true,
-            "type": "boolean"
+            "anyOf": [
+                {
+                    "type": "null"
+                },
+                {
+                    "type": "number"
+                }
+            ],
+            "default": null
         }
     },
     "required": [
         "algorithm",
         "saveat",
         "adaptive",
-        "dt",
+        "force_dtmin",
         "abstol",
         "reltol",
         "maxiters",

python/ribasim/ribasim/config.py

@@ -38,15 +38,18 @@ class DiscreteControl(BaseModel):
 
 
 class Solver(BaseModel):
-    reltol: float = 0.001
-    saveat: Union[List[float], float] = []
-    maxiters: int = 1000000000
     autodiff: bool = True
     adaptive: bool = True
+    force_dtmin: bool = False
+    dt: Optional[float] = None
+    saveat: Union[List[float], float] = []
+    maxiters: int = 1000000000
+    dtmin: Optional[float] = None
+    sparse: bool = True
     algorithm: str = "QNDF"
     abstol: float = 1e-06
-    dt: float = 0
-    sparse: bool = True
+    reltol: float = 0.001
+    dtmax: Optional[float] = None
 
 
 class FlowBoundary(BaseModel):
@@ -131,7 +134,10 @@ class Config(BaseModel):
                 "algorithm": "QNDF",
                 "saveat": [],
                 "adaptive": True,
-                "dt": 0,
+                "dt": None,
+                "dtmin": None,
+                "dtmax": None,
+                "force_dtmin": False,
                 "abstol": 1e-06,
                 "reltol": 0.001,
                 "maxiters": 1000000000,