Numerical considerations

We want to solve the following initial value problem: \[ \begin{cases} \frac{\text{d}\mathbf{u}}{\text{d}t} = \mathbf{f}(\mathbf{u},t) \quad t_0 < t < t_\text{end} \\ \mathbf{u}(t_0) = \mathbf{u}_0 \end{cases}, \tag{1}\]

where \(\mathbf{f}\) denotes water_balance! and \(\mathbf{u_0} = \mathbf{0}\) the initial cumulative flows (and the PID integrals which also start out at \(0\)).

In general \(\mathbf{f}\) is a non-linear function in \(\mathbf{u}\). These non-linearities are introduced by e.g.:

The problem Equation 1 can be solved by various numerical time-integration methods. To do this the time interval \([t_0,t_\text{end}]\) is discretized into a finite number of time points \(t_0 < t_1 < \ldots < t_N = t_\text{end}\) for which approximate solutions \(\mathbf{w}_n \approx \mathbf{u}(t_n)\) are computed. In general we do not assume a fixed timestep (the interval between successive points in time). Rather, the solver attempts to make as large a step as possible while keeping error tolerances within requirements. The solver settings section details the available configuration options.

1 Example numerical methods

This section discusses two relatively simple numerical methods, mainly to demonstrate the difference between explicit and implicit methods.

1.1 Euler forward

The simplest numerical method is Euler forward: \[ \mathbf{w}_{n+1} = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_n, t_n). \tag{2}\]

Here \(\mathbf{w}_{n+1}\) is given as a simple explicit function of \(\mathbf{w}_n\).

1.2 Euler backward

Euler backward is formulated as follows: \[ \mathbf{w}_{n+1} = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1},t_{n+1}). \tag{3}\]

Note that this is an implicit equation for \(\mathbf{w}_{n+1}\), which is non-linear because of the non-linearity of \(\mathbf{f}\).

Generally one of the following iterative methods is used for finding solutions to non-linear equations like this:

  • Picard iteration for fixed points. This method aims to approximate \(\mathbf{w}_{n+1}\) as a fixed point of the function \[ \mathbf{g}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) \] by iterating \(\mathbf{g}\) on an initial guess of \(\mathbf{w}_{n+1}\);
  • Newton iterations: approximate \(\mathbf{w}_{n+1}\) as a root of the function \[ \mathbf{h}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) - \mathbf{x}, \] by iteratively finding the root of its linearized form:

\[\begin{align} \mathbf{0} =& \mathbf{h}(\mathbf{w}_{n+1}^k) + \mathbf{J}(\mathbf{h})(\mathbf{w}_{n+1}^k)(\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k) \\ =& \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1}^k,t_{n+1}) - \mathbf{w}_{n+1}^k \\ +&\left[(t_{n+1}-t_n)\mathbf{J}(\mathbf{f})(\mathbf{w}_{n+1}^k)-\mathbf{I}\right](\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k). \end{align}\] Note that this thus requires an evaluation of the Jacobian of \(\mathbf{f}\) and solving a linear system per iteration.

2 The advantage of implicit methods

The implicit method Equation 3 is a coupled system of equations for \(\mathbf{w}_{n+1}\), while the explicit Equation 2 is fully decoupled. This means in general that in a time integration step with an implicit method the basins communicate information with eachother, while in an explicit method this does not happen. A consequence of this is that local events (e.g. a pump turns on) propagate slowly through the model using an explicit method but quickly using an implicit method, making implicit methods more stable.

3 The Jacobian

The Newton iterations above require the Jacobian of \(\mathbf{f}\), an \(N \times N\) matrix where \(N\) is the number of states. Its entries are given by \[ J_{i,j} = \frac{\partial f_i}{\partial u_j}, \] i.e. \(J_{i,j}\) quantifies how \(f_i\), the time derivative of state \(i\), changes with respect to changes in state \(j\). Most of these entries are \(0\), because flows in distant parts of the model do not depend on each other.

There are several ways to compute the Jacobian:

  • Using finite difference methods; however these are relatively expensive to compute and introduce truncation errors.
  • Hardcoding analytical expressions for the partial derivatives; these have the advantages over finite difference methods that they are cheaper to compute and don’t introduce truncation errors. However, computing and coding these partial derivatives is complex and thus error prone, and requires maintenance for every adjustment of \(\mathbf{f}\).
  • Using automatic differentiation; this has the advantages of being fast to compute and not requiring much extra implementation other than implementing \(\mathbf{f}\) in a way that is compatible with an automatic differentiation package, in the case of Ribasim that is ForwardDiff.jl.

Ribasim uses automatic differentiation by default. With forward mode automatic differentiation, one evaluation of \(\mathbf{f}\) yields one column of the Jacobian, or more precisely the Jacobian applied to one chosen direction. A naive computation of the Jacobian thus costs \(N\) evaluations of \(\mathbf{f}\), which is prohibitively expensive for large models.

3.1 Sparse matrix coloring

If two columns of the Jacobian have no nonzero entries in the same row, they can be computed in a single evaluation of \(\mathbf{f}\) by perturbing both corresponding states at once: the contributions can be told apart afterwards based on the sparsity pattern. The columns of a sparse Jacobian can be partitioned into groups (colors) of such structurally orthogonal columns, and only one evaluation per color is needed. Ribasim detects the sparsity pattern automatically with SparseConnectivityTracer.jl and computes the coloring with SparseMatrixColorings.jl.

The number of colors depends on the structure of the Jacobian. For the Jacobian of \(\mathbf{f}\) this structure is unfavorable: a perturbation of a cumulative flow state changes the storage of the Basins it connects, and that affects all flows connected to these Basins. Hence all flow states connected to the same Basin end up in different colors, and the number of colors is at least the number of connections of the most connected Basin. Therefore Ribasim does not compute the Jacobian of \(\mathbf{f}\) directly, but the Jacobian of an intermediate function which can be colored with very few colors, as explained below.

3.2 Structure of the right hand side

We split the state vector into the cumulative flows \(\mathbf{u}^\text{flow}\) and the PID integrals \(\mathbf{u}^\text{pid}\). The Basin storages depend linearly on the cumulative flows (see the equations): \[ \mathbf{S} = \mathbf{S}_0(t) + M\mathbf{u}^\text{flow}, \] where \(M\) is the \(n_\text{basin} \times n_\text{flow}\) signed incidence matrix: \(M_{b,k} = -1\) if Basin \(b\) is upstream of flow \(k\), \(M_{b,k} = 1\) if Basin \(b\) is downstream of flow \(k\) and \(M_{b,k} = 0\) otherwise. Evaporation and infiltration are flows with their Basin upstream and nothing downstream.

For every flow \(k\) we then take a copy of the storage of its upstream Basin \(S^\text{up}_k\) and its downstream Basin \(S^\text{down}_k\): \[ \mathbf{S}^\text{up} = P^\text{up}\mathbf{S}, \quad \mathbf{S}^\text{down} = P^\text{down}\mathbf{S}, \] where \(P^\text{up}\) and \(P^\text{down}\) are constant selection matrices (note that \(M = (P^\text{down} - P^\text{up})^T\)). The flows are computed from these storage copies, the PID integrals \(\mathbf{u}^\text{pid}\) and the compound variables \(\mathbf{c}\) of the ContinuousControl nodes: \[ \mathbf{f}^\text{flow} = \mathbf{q}(\mathbf{v}), \quad \mathbf{v} = \left(\mathbf{S}^\text{up}, \mathbf{S}^\text{down}, \mathbf{u}^\text{pid}, \mathbf{c}\right). \]

The Jacobian \(\partial\mathbf{q}/\partial\mathbf{v}\) is what Ribasim computes with automatic differentiation. It is not the Jacobian of \(\mathbf{f}\), but it is very cheap to compute:

  • \(\partial\mathbf{q}/\partial\mathbf{S}^\text{up}\) and \(\partial\mathbf{q}/\partial\mathbf{S}^\text{down}\) are diagonal, since flow \(k\) only depends on its own storage copies \(S^\text{up}_k\) and \(S^\text{down}_k\). The exception are flows controlled by a PidControl node with a derivative term, which depend on all flows of the listened Basin;
  • \(\partial\mathbf{q}/\partial\mathbf{u}^\text{pid}\) and \(\partial\mathbf{q}/\partial\mathbf{c}\) have at most one nonzero per row, for the flows that are PID controlled or continuously controlled respectively.

These blocks can be colored with only a handful of colors, independent of the size of the model. The Jacobian of \(\mathbf{f}\) follows from the chain rule, where only the multiplications with the constant matrices \(P^\text{up}M\) and \(P^\text{down}M\) remain: \[ \frac{\partial \mathbf{f}^\text{flow}}{\partial \mathbf{u}^\text{flow}} = Q_S M, \quad Q_S = \frac{\partial\mathbf{q}}{\partial\mathbf{S}^\text{up}}P^\text{up} + \frac{\partial\mathbf{q}}{\partial\mathbf{S}^\text{down}}P^\text{down}, \] (ignoring control for now), where \(Q_S = \partial\mathbf{q}/\partial\mathbf{S}\) is the \(n_\text{flow} \times n_\text{basin}\) matrix with the dependence of the flows on the Basin storages. This matrix is never assembled; instead the structure is exploited directly in the linear solve.

3.3 The linear solve

Every Newton iteration of an implicit method requires solving a linear system of the form \[ \left(J - \gamma^{-1}I\right)\mathbf{x} = \mathbf{b}, \] where \(\gamma\) is proportional to the timestep (compare the Euler backward iteration above, where \(\gamma = t_{n+1} - t_n\)). This system has a dimension equal to the number of states, which is much larger than the number of Basins. Using \(J = Q_S M\) and defining \(\mathbf{s} = M\mathbf{x}\), the system reads \[ -\gamma^{-1}\mathbf{x} + Q_S\mathbf{s} = \mathbf{b} \quad \Rightarrow \quad \mathbf{x} = \gamma\left(Q_S\mathbf{s} - \mathbf{b}\right). \tag{4}\] Multiplying the left equation by \(M\) gives a system in terms of \(\mathbf{s}\) only: \[ \left(J_\text{inner} - \gamma^{-1}I_{n_\text{basin}}\right)\mathbf{s} = M\mathbf{b}, \quad J_\text{inner} = M Q_S. \tag{5}\] Hence the linear solve is done in two steps:

  1. Solve the \(n_\text{basin} \times n_\text{basin}\) system Equation 5 for the storage changes \(\mathbf{s}\). The sparsity pattern of \(J_\text{inner}\) is that of the adjacency matrix of the Basins, which is very sparse; by default this system is solved with a sparse LU factorization.
  2. Compute the solution in terms of the states \(\mathbf{x}\) with Equation 4, which is just a cheap matrix-vector product.

3.3.1 PID control

The PID integral state of a PidControl node \(j\) evolves with the error \(f^\text{pid}_j = e_j = h^\text{target}_j - h(S_{l_j})\), where \(l_j\) is the listened Basin, so \[ \frac{\partial f^\text{pid}_j}{\partial \mathbf{u}^\text{flow}} = -\frac{1}{A_{l_j}}M_{l_j,:}, \] where \(A_{l_j}\) is the area of the listened Basin. The listened Basin is directly up- or downstream of the controlled Pump or Outlet, so the dependence of the controlled flow on the listened Basin storage through the proportional term is already included in \(Q_S\). Let \(G = \partial\mathbf{q}/\partial\mathbf{u}^\text{pid}\) and let \(E\) be the matrix that selects the listened Basin per PidControl node and divides by its area. The PID rows of the linear system then give \[ \mathbf{x}^\text{pid} = -\gamma\left(\mathbf{b}^\text{pid} + E\mathbf{s}\right), \] and substituting this into the flow rows leads to \[ \mathbf{x}^\text{flow} = \gamma\left(Q_S\mathbf{s} + G\mathbf{x}^\text{pid} - \mathbf{b}^\text{flow}\right), \] with the inner system \[ \left(M Q_S - \gamma M G E - \gamma^{-1}I_{n_\text{basin}}\right)\mathbf{s} = M\mathbf{b}^\text{flow} + \gamma M G \mathbf{b}^\text{pid}. \] So the PID integral states are eliminated in the same way as the flow states, and the dimension of the inner system remains the number of Basins.

3.3.2 Continuous control

The compound variable \(c_m\) of a ContinuousControl node depends on Basin storages and on (listened) flows. Their derivatives \(\partial\mathbf{c}/\partial\mathbf{S}\) and \(\partial\mathbf{c}/\partial\mathbf{q}\) are also computed with automatic differentiation. For a continuously controlled flow this adds the term \[ \frac{\partial\mathbf{q}}{\partial\mathbf{c}}\left(\frac{\partial\mathbf{c}}{\partial\mathbf{S}} + \frac{\partial\mathbf{c}}{\partial\mathbf{q}}Q_S\right) \] to \(Q_S\), and thus \(M\) times this term to \(J_\text{inner}\). Here only the direct dependence of listened flows on their own Basin storages is taken into account, since chaining control (listening to a flow that is itself controlled by PidControl or ContinuousControl) is not supported.

3.3.3 Exactness

The linear solve above is equivalent to solving the full linear system with the Jacobian of \(\mathbf{f}\), as long as control is not chained. If it is, the Jacobian used in the linear solve is approximate. An approximate Jacobian only affects the convergence speed of the Newton iterations, not the accuracy of the converged solution, since the latter is determined by \(\mathbf{f}\) itself.

The reduced linear solve can be disabled with the solver setting reduced_implicit_solve = false, see the solver settings. In that case the (sparse) Jacobian of \(\mathbf{f}\) is computed directly with automatic differentiation and sparse matrix coloring, and the full linear system is solved. This is mainly useful for verification, as it is considerably slower for larger models.

4 Continuity considerations

The convergence of the Newton method can be proven given certain properties of \(\mathbf{f}\) around the initial guess and root to find. An important aspect is the smoothness of \(\mathbf{f}\). The basin profiles and \(Q(h)\) relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.

Control mechanisms can change parameters of \(\mathbf{f}\) discontinuously, leading to discontinuities of \(\mathbf{f}\). This however does not yield problems for the time integration methods in DifferentialEquations.jl, since the callback mechanisms used to change these parameters make the solver take these discontinuities into account.

4.1 Basin profiles

The basin profiles affect \(\mathbf{f}\) in many ways, anywhere where a basin level or area is required.

Note

This section needs to be updated and extended after once this issue is resolved.

4.2 Q(h) relations

TabulatedRatingCurve nodes contribute to \(\mathbf{f}\) with terms of the following form:

\[ Q(h(u)) \]

where the continuity of this term is given by the least continuous of \(Q\) and \(h\).

4.3 Empty basins

Reduction factors are introduced at several points in the definition of \(\mathbf{f}\) to smooth out otherwise discontinuous transitions (e.g. the flow rate of a pump going to zero when the source basin dries out). If flows are not too large with respect to basin storage, this will prevent basins from reaching 0. Rather, the basin gets a very small storage. The reduction factors help with performance, but are also an important tool to avoid getting negative storage in basins. Negative storage needs to be avoided since it is not a real solution, and would introduce water into the model that doesn’t exist. Another tool used to avoid negative storage is the isoutoutofdomain option, which Ribasim makes use of. This rejects timesteps that lead to negative storage, instead retrying with a smaller timestep.

5 Error control for cumulative flows

The states in \(\mathbf{u}\) are cumulative quantities that grow monotonically over time. With the default SciML error control the tolerance per state would be roughly \(\text{abstol} + |\mathbf{u}| * \text{reltol}\), see Basic Stepsize Control. That means the tolerance is based on the cumulative flow so far rather than on the flow over the current timestep, so it keeps loosening as the simulation progresses: after a thousand seconds at \(1\text{ m}^3/\text{s}\) an effective reltol of \(10^{-5}\) has degraded to \(10^{-2}\). It also makes the tolerance depend on arbitrary offsets: a large Basin with little throughflow would get a very loose tolerance purely because of its storage.

Ribasim therefore overrides calculate_residuals so that each state is scaled by the magnitude of its change over the timestep instead of its absolute magnitude:

\[ \text{residual}_i = \frac{\tilde{u}_i}{\text{abstol} + |u_{1,i} - u_{0,i}| * \text{reltol}} \]

This keeps the meaning of abstol and reltol constant over the whole simulation, and makes them independent of the size of the Basin storages. The same scaling is applied to the increment of the nonlinear solver, so that the nonlinear solve and the error estimate weigh the states in the same way.

Note that this changes the units in which the tolerances should be interpreted. abstol is an absolute tolerance on the volume passed by a state during a single timestep (in \(\text{m}^3\)), and reltol is relative to that same per-step volume. It acts as a floor for states that carry almost no flow during a timestep, where the reltol term vanishes; setting it too low makes such states dictate the timestep for the whole model. Since it is a volume per timestep and not a rate, abstol should not be scaled with dt.

The default of \(10^{-4}\ \text{m}^3\) is the largest value that is accurate across model scales, from a handful of Basins up to a national model. Because it is an absolute volume it does not scale with the model, so a large model can often raise it to \(10^{-2}\ \text{m}^3\) for a substantial speedup at no measurable loss of accuracy, while the same setting on a small model is roughly two orders of magnitude less accurate. If you raise it, verify against a run with a tighter tolerance.