From 94c79e02b787dd3178d93905f5ab4f2fec894b7f Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Thu, 5 Oct 2023 11:55:38 +0200 Subject: [PATCH 1/2] start using isoutofdomain (#639) This rejects any steps that lead to a negative storage. See also: https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#My-ODE-goes-negative-but-should-stay-positive,-what-tools-can-help? When using adaptive solvers this will cause the dt to go down. If it is not adaptive and implicit, it will fail: ``` Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt. ``` Both options seem desired behavior. I managed to still get negative storages with an explicit solver with a large fixed dt, and a flux that empties a basin in well under a timestep. But adaptive is the default and should generally be used. The `du` correction in `water_balance!` is also removed as it should normally not be needed anymore. And if it does, I'd prefer to not distort the water balance. One alternative is to use the [`PositiveDomain`](https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain) callback, but this only works for `adaptive=false`, and is not easily adapted to enforce only the storage part of the state to be positive. Also adds an unrelated test to avoid saving a t more than once when not needed. --- core/src/bmi.jl | 1 + core/src/solve.jl | 8 -------- core/test/run_models.jl | 1 + 3 files changed, 2 insertions(+), 8 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index fdfda2ef5..3eb0278ea 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -126,6 +126,7 @@ 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, diff --git a/core/src/solve.jl b/core/src/solve.jl index fa62f3582..6b83d54f2 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -1154,14 +1154,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 diff --git a/core/test/run_models.jl b/core/test/run_models.jl index 4f0b37301..6a59c0cc6 100644 --- a/core/test/run_models.jl +++ b/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 From df2099fc627422fad6f53e0d7afbe229d320bf10 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Thu, 5 Oct 2023 14:09:48 +0200 Subject: [PATCH 2/2] add more solver options for timestepping (#641) https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/ This provides options for bounding the allowed dt for adaptive timestepping, even forcing it to continue despite dtmin being too large, which obviously should be used with caution. --- core/src/bmi.jl | 12 +++++- core/src/config.jl | 16 +++++--- core/test/docs.toml | 21 ++++++----- docs/core/usage.qmd | 6 ++- docs/schema/Config.schema.json | 5 ++- docs/schema/Solver.schema.json | 64 +++++++++++++++++++++++++------- python/ribasim/ribasim/config.py | 18 ++++++--- 7 files changed, 104 insertions(+), 38 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 3eb0278ea..c35e820ed 100644 --- a/core/src/bmi.jl +++ b/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; @@ -129,7 +134,10 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model 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, diff --git a/core/src/config.jl b/core/src/config.jl index a677a33b9..376e3ffa5 100644 --- a/core/src/config.jl +++ b/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 diff --git a/core/test/docs.toml b/core/test/docs.toml index 67fc33019..b3335c1fa 100644 --- a/core/test/docs.toml +++ b/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 diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 54902ba12..c66445f7b 100644 --- a/docs/core/usage.qmd +++ b/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. diff --git a/docs/schema/Config.schema.json b/docs/schema/Config.schema.json index e35529df9..290954732 100644 --- a/docs/schema/Config.schema.json +++ b/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, diff --git a/docs/schema/Solver.schema.json b/docs/schema/Solver.schema.json index ca6c57bd9..fc95a4177 100644 --- a/docs/schema/Solver.schema.json +++ b/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", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 185f98877..3ee4ff0a9 100644 --- a/python/ribasim/ribasim/config.py +++ b/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,