Debugging models

1 Unexpected flow rates

When a connector node (e.g., TabulatedRatingCurve, Pump, Outlet, ManningResistance) produces unexpected flow rates, it’s important to understand how Ribasim calculates flows and which reduction factors may apply.

1.1 Review the equations

Each connector node type has specific equations that govern its behavior. Review the equations for your node type in the reference documentation. For example, see the TabulatedRatingCurve equations.

1.2 Check reduction factors

Ribasim applies several reduction factors that can limit flows. To verify that none of them unexpectedly applies, check the following examples:

Upstream Basin level and storage: If the upstream Basin is nearly empty, a reduction factor prevents extracting more water than available.
Downstream Basin level: If the downstream Basin level is higher, only Pumps can flow against gravity.
Connector node parameters: Review the node’s parameters in your model input to ensure they match your expectations. Parameters like max_flow_rate or min_upstream_level directly affect the computed flow.

1.3 Check control and allocation

If the node is subject to control or allocation:

Discrete Control: Check the control results to see which control state is active and when it changes. The active state determines which parameter set is used.
Allocation: Check the allocation results to see the allocated flow for the node.

In all these cases the reduction factors still apply, and may be the reason allocated flows cannot be supplied.

2 Slow models

When your model is slow, it is usually only a handful of nodes that are hard to solve. Ribasim points them out in three places, from coarse to detailed.

2.1 1. Read the convergence bottlenecks in the log

At the end of every run, whether it succeeded, failed or was interrupted, Ribasim lists the nodes that were hardest to solve:

┌ Info: Convergence bottlenecks in descending order of severity (1.0 is the worst):
│   Outlet #2310 = 1.0
│   Pump #10720 = 0.65
│   ManningResistance #2383 = 0.61
│   ManningResistance #5494 = 0.36
└   TabulatedRatingCurve #7185 = 0.36

The value is a severity score: the share of the local error estimate that Ribasim attributes to that node, averaged over time and normalized so that the worst node of every timestep scores 1. It ranks the nodes, it does not measure how wrong they are. A high value does not mean the node is misconfigured or that its results are inaccurate; it means the node is the one holding the solver back. The top entry is therefore always at or near 1, and it is the node to look at first. Because the solver uses one timestep for the whole model, a single node that is hard to solve slows down everything else with it. Unphysical parameters are a common cause: in the example above, Outlet #2310 had a flow_rate of 99999.0, and ManningResistance #2383 a roughness low enough to hardly have any friction.

2.2 2. Find out when and where, using the results

The same measure is written per node and per saveat interval as the convergence variable in basin.nc and flow.nc, see the output files. Plot it over time to see whether a node is always problematic or only during, for example, a dry spell or a control switch.

2.3 3. Find out why, using the solver statistics

solver_stats.nc records per saveat interval how much work the solver did and why it had to redo it:

variable	meaning
`dt`	average calculation timestep size in seconds; this is the number to grow
`accepted_timesteps`	timesteps that were kept
`rejected_timesteps`	timesteps that were thrown away and retried with a smaller step
`rejected_local_error`	rejected because the estimated error over the step was too large
`rejected_nonlinear_solve`	rejected because the solver could not solve the equations to the required accuracy
`rejected_out_of_domain`	rejected because a Basin would have run dry within the step
`mean_order`	average order of the accepted timesteps
`rhs_calls`, `linear_solves`	raw solver work
`computation_time`	wall time in milliseconds spent on the interval

A high rejected_local_error means the model is stiff: some node reacts far faster than the timestep you would like to take. Combine it with the convergence bottlenecks above to find which node, then look at its parametrization.

A high rejected_nonlinear_solve or rejected_out_of_domain instead points at a Basin that is (nearly) empty, or at a discontinuity such as a control switch.

A mean_order that is stuck around 1 or 2 while the solver can go up to 5 is a second sign of stiffness or of a discontinuity that keeps resetting the solver.

2.4 What to change

Once you know the node, the usual causes are:

Unphysical parameter values, often placeholders such as 9999 for flow_rate. These are not harmless: the flow rate multiplies the slope of the flow-versus-level relation, and that slope is what sets the timestep.
A small Basin bottom area, so even a modest flow changes its level quickly. Start its profile with an area of at least 10 m².
A ManningResistance between two Basins at almost the same level, which sits on the steep part of its equation. Check its profile_width, length and manning_n, or merge the two Basins.

Loosening abstol and reltol in the solver settings can speed up a model at the cost of accuracy. Compare the results with a run using stricter tolerances to check whether the differences are acceptable.

3 Unstable models

When your model exits with a message like so:

┌ Error: The model exited at model time 2024-01-27T14:46:17.791 with return code Unstable. See https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#retcodes

the solver kept shrinking the timestep until it gave up. The convergence bottlenecks are logged in this case as well, and they name the node to start with.

To see what happened just before the model gave up, rerun with saveat = 0 in the solver settings. All output up to the moment of the error is written, so the last few rows of basin.nc, flow.nc and solver_stats.nc show the state of the model as it went unstable. The model might then instead exit with

┌ Error: Too large water balance error
│   id = Basin #2
│   balance_error = 0.0017985748886167501
│   relative_error = 1.3503344464431657

which helps you pin down the problematic node(s).

The rejected_* variables in solver_stats.nc tell you which wall the solver ran into:

rejected_out_of_domain rising means a Basin is being drained faster than it can supply. Check the outgoing connector nodes and the Basin profile.
rejected_nonlinear_solve rising means the equations themselves became hard to solve, typically around a discontinuity such as a control switch or a TabulatedRatingCurve with a sharp kink.
rejected_local_error rising means the model is stiff rather than truly unstable; see slow models.