Introduce config option for steady state at simulation start

core/src/config.jl

@@ -92,6 +92,7 @@ const nodetypes = collect(keys(nodekinds))
 
 @option struct Solver <: TableOption
     algorithm::String = "QNDF"
+    steady_start::Bool = true
     saveat::Float64 = 86400.0
     dt::Union{Float64, Nothing} = nothing
     dtmin::Float64 = 0.0

core/src/model.jl

@@ -129,10 +129,15 @@ function Model(config::Config)::Model
     jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing
     RHS = ODEFunction(water_balance!; jac_prototype)
 
-    @timeit_debug to "Compute steady start state" begin
-        steady_state_start!()
+    if config.solver.steady_start
+        @timeit_debug to "Compute steady start state" begin
+            # Look for steady state on first simulated day
+            sol =
+                solve(SteadyStateProblem(RHS, u0, parameters), DynamicSS(; tspan = 86400.0))
+            u0 .= sol.u
+        end
+        @debug "Compute steady start state"
     end
-    @debug "Compute steady start state"
 
     @timeit_debug to "Setup ODEProblem" begin
         prob = ODEProblem(RHS, u0, timespan, parameters)

core/src/util.jl

@@ -814,10 +814,3 @@ function set_previous_flows!(integrator)
     copyto!(flow_prev, flow)
     copyto!(vertical_flux_prev, vertical_flux)
 end
-
-function steady_state_start(ode_prob)
-    # Find steady state in first day
-    sol = solve(SteadyStateProblem(ode_prob), DynamicSS(; tspan = 86400.0))
-    display(propertynames(sol))
-    sol
-end

core/test/docs.toml

@@ -26,6 +26,7 @@ use_allocation = false             # optional, default false
 
 [solver]
 algorithm = "QNDF"  # optional, default "QNDF"
+steady_start = true # optional, default true
 saveat = 86400      # optional, default 86400, 0 saves every timestep, inf saves only at start- and endtime
 dt = 60.0           # optional, remove for adaptive time stepping
 dtmin = 0.0         # optional, default 0.0

docs/reference/usage.qmd

@@ -30,6 +30,8 @@ With adaptive timestepping, `dtmin` and `dtmax` control the minimum and maximum
 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.
 
+An initial state is provided in the `Basin / state` table. If `steady_start = true` then the solver first searches for a steady state solution at the start of the simulation (with the provided state as initial condition).
+
 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`.