From f018c038474becbcf4d22ef90239b5046ee83ab0 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Thu, 3 Oct 2024 17:56:40 +0200 Subject: [PATCH 1/3] Add more solvers --- core/src/config.jl | 14 +++++++++----- core/test/run_models_test.jl | 2 +- docs/reference/usage.qmd | 2 +- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/core/src/config.jl b/core/src/config.jl index 1aa849892..2b31d3939 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -18,8 +18,8 @@ using OrdinaryDiffEqNonlinearSolve: NLNewton using OrdinaryDiffEqLowOrderRK: Euler, RK4 using OrdinaryDiffEqTsit5: Tsit5 using OrdinaryDiffEqSDIRK: ImplicitEuler, KenCarp4, TRBDF2 -using OrdinaryDiffEqBDF: QNDF -using OrdinaryDiffEqRosenbrock: Rodas5, Rosenbrock23 +using OrdinaryDiffEqBDF: FBDF, QNDF +using OrdinaryDiffEqRosenbrock: Rosenbrock23, Rodas4P, Rodas5P export Config, Solver, Results, Logging, Toml export algorithm, @@ -213,20 +213,24 @@ Map from config string to a supported algorithm type from [OrdinaryDiffEq](https Supported algorithms: - `QNDF` +- `FBDF` - `Rosenbrock23` - `TRBDF2` -- `Rodas5` +- `Rodas4P` +- `Rodas5P` - `KenCarp4` - `Tsit5` - `RK4` - `ImplicitEuler` - `Euler` """ -const algorithms = Dict{String, Type}( +algorithms = Dict{String, Type}( "QNDF" => QNDF, + "FBDF" => FBDF, "Rosenbrock23" => Rosenbrock23, "TRBDF2" => TRBDF2, - "Rodas5" => Rodas5, + "Rodas4P" => Rodas4P, + "Rodas5P" => Rodas5P, "KenCarp4" => KenCarp4, "Tsit5" => Tsit5, "RK4" => RK4, diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 065df785b..20cc0b492 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -280,7 +280,7 @@ end @test sparse_fdm.integrator.u ≈ sparse_ad.integrator.u atol = 4 @test dense_fdm.integrator.u ≈ sparse_ad.integrator.u atol = 4 - config = Ribasim.Config(toml_path; solver_algorithm = "Rodas5", solver_autodiff = true) + config = Ribasim.Config(toml_path; solver_algorithm = "Rodas5P", solver_autodiff = true) time_ad = Ribasim.run(config) @test successful_retcode(time_ad) @test time_ad.integrator.u ≈ sparse_ad.integrator.u atol = 10 diff --git a/docs/reference/usage.qmd b/docs/reference/usage.qmd index 354a3afc0..50f1f63cd 100644 --- a/docs/reference/usage.qmd +++ b/docs/reference/usage.qmd @@ -21,7 +21,7 @@ If your model does not converge, or your performance is lower than expected, it 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 timestepping. -The full list of available solvers is: `QNDF`, `Rosenbrock23`, `TRBDF2`, `Rodas5`, `KenCarp4`, `Tsit5`, `RK4`, `ImplicitEuler`, `Euler`. +The full list of available solvers is: `QNDF`, `FBDF`, `Rosenbrock23`, `Rodas4P`, `Rodas5P`, `TRBDF2`, `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/). By default Ribasim uses adaptive timestepping, though not all algorithms support adaptive timestepping. From 0c0a00c0859690e6f78298098b26b572c2302ae9 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Thu, 3 Oct 2024 18:02:33 +0200 Subject: [PATCH 2/3] Keep it const --- core/src/config.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/config.jl b/core/src/config.jl index 2b31d3939..0ae043004 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -224,7 +224,7 @@ Supported algorithms: - `ImplicitEuler` - `Euler` """ -algorithms = Dict{String, Type}( +const algorithms = Dict{String, Type}( "QNDF" => QNDF, "FBDF" => FBDF, "Rosenbrock23" => Rosenbrock23, From 11aa5c35addb8cee5e2da762625957014084a6a8 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Thu, 3 Oct 2024 18:21:45 +0200 Subject: [PATCH 3/3] Try Rodas5P default solver --- core/src/config.jl | 2 +- core/test/config_test.jl | 2 +- core/test/docs.toml | 2 +- docs/concept/equations.qmd | 2 +- docs/reference/usage.qmd | 2 +- docs/references.bib | 11 ----------- python/ribasim/ribasim/config.py | 4 ++-- python/ribasim/tests/test_model.py | 2 +- 8 files changed, 8 insertions(+), 19 deletions(-) diff --git a/core/src/config.jl b/core/src/config.jl index 0ae043004..67a5548b7 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -94,7 +94,7 @@ end const nodetypes = collect(keys(nodekinds)) @option struct Solver <: TableOption - algorithm::String = "QNDF" + algorithm::String = "Rodas5P" saveat::Float64 = 86400.0 dt::Union{Float64, Nothing} = nothing dtmin::Float64 = 0.0 diff --git a/core/test/config_test.jl b/core/test/config_test.jl index 8f9cb3e8f..87e070c7b 100644 --- a/core/test/config_test.jl +++ b/core/test/config_test.jl @@ -41,7 +41,7 @@ end using Ribasim: convert_saveat, convert_dt, Solver, algorithm solver = Solver() - @test solver.algorithm == "QNDF" + @test solver.algorithm == "Rodas5P" Solver(; algorithm = "Rosenbrock23", autodiff = true, diff --git a/core/test/docs.toml b/core/test/docs.toml index 39b7b2894..0a9db1657 100644 --- a/core/test/docs.toml +++ b/core/test/docs.toml @@ -25,7 +25,7 @@ timestep = 86400 # optional (required if use_allocation = true use_allocation = false # optional, default false [solver] -algorithm = "QNDF" # optional, default "QNDF" +algorithm = "Rodas5P" # optional, default "Rodas5P" 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 diff --git a/docs/concept/equations.qmd b/docs/concept/equations.qmd index fdd88df4d..044a0cdcc 100644 --- a/docs/concept/equations.qmd +++ b/docs/concept/equations.qmd @@ -160,7 +160,7 @@ There are many things that can influence the calculations times, for instance: - [Solver tolerance](https://diffeq.sciml.ai/stable/basics/faq/#What-does-tolerance-mean-and-how-much-error-should-I-expect): By default we use absolute and relative tolerances of `1e-6` and `1e-5` respectively. - [ODE solvers](https://diffeq.sciml.ai/stable/solvers/ode_solve/): - The `QNDF` method we use is robust to oscillations and massive stiffness, however other solvers should be tried as well. + The `Rodas5P` method we use is robust to oscillations and massive stiffness, however other solvers should be tried as well. - Forcing: Every time new forcing data is injected into the model, it needs to pause. Moreover, the larger the forcing fluxes are, the bigger the shock to the system, leading to smaller timesteps and thus longer simulation times. diff --git a/docs/reference/usage.qmd b/docs/reference/usage.qmd index 50f1f63cd..df7f75f9a 100644 --- a/docs/reference/usage.qmd +++ b/docs/reference/usage.qmd @@ -19,7 +19,7 @@ The solver section in the configuration file is entirely optional, since we aim 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]. +The default solver is `algorithm = "Rodas5P"`, which is a 5th order A-stable stiffly stable Rosenbrock method. It is an implicit method that supports the default adaptive timestepping. The full list of available solvers is: `QNDF`, `FBDF`, `Rosenbrock23`, `Rodas4P`, `Rodas5P`, `TRBDF2`, `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/). diff --git a/docs/references.bib b/docs/references.bib index 2fd117c89..00efd9330 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -87,14 +87,3 @@ @misc {pdoktopnl url = "https://www.pdok.nl/downloads/-/article/basisregistratie-topografie-brt-topnl", note = "[Online; accessed 31-August-2022]" } - -@article{shampine1997matlab, - title={The matlab ode suite}, - author={Shampine, Lawrence F and Reichelt, Mark W}, - journal={SIAM journal on scientific computing}, - volume={18}, - number={1}, - pages={1--22}, - year={1997}, - publisher={SIAM} -} diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 2511c1eac..ccb415ce9 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -86,7 +86,7 @@ class Solver(ChildModel): Attributes ---------- algorithm : str - The used numerical time integration algorithm (Optional, defaults to QNDF) + The used numerical time integration algorithm (Optional, defaults to Rodas5P) saveat : float Time interval in seconds between saves of output data. 0 saves every timestep, inf only saves at start- and endtime. (Optional, defaults to 86400) @@ -111,7 +111,7 @@ class Solver(ChildModel): Whether automatic differentiation instead of fine difference is used to compute the Jacobian. (Optional, defaults to true) """ - algorithm: str = "QNDF" + algorithm: str = "Rodas5P" saveat: float = 86400.0 dt: float | None = None dtmin: float | None = None diff --git a/python/ribasim/tests/test_model.py b/python/ribasim/tests/test_model.py index ba703d5fa..86e88776d 100644 --- a/python/ribasim/tests/test_model.py +++ b/python/ribasim/tests/test_model.py @@ -29,7 +29,7 @@ def test_repr(basic): def test_solver(): solver = Solver() - assert solver.algorithm == "QNDF" # default + assert solver.algorithm == "Rodas5P" # default assert solver.saveat == 86400.0 solver = Solver(saveat=3600.0)