Try Rodas5P as default solver

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

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,

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

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.
 

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/).

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}
-}

python/ribasim/ribasim/config.py

@@ -85,7 +85,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)
@@ -110,7 +110,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

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)