This is a bit out of date with #3113, but the general idea is still valid.
Let's investigate the structure of the linear system that emerges in #3113 for $x_s$. The matrix there can be written as:
$$
I_n - \gamma \left[S_\text{down}^T - S_\text{up}^T\right]\left[F + \frac{\partial q}{\partial c}\frac{\partial c}{s}\right],
$$
where
$$
\begin{align}
F &=& \text{diag}(v_1)S_\text{up}\text{diag}(w_1) + \\
&&\text{diag}(v_2)S_\text{up}\text{diag}(w_2) + \\
&&\text{diag}(v_3)S_\text{down}\text{diag}(w_1) + \\
&&\text{diag}(v_4)S_\text{down}\text{diag}(w_2) \\
S_\text{PID} &=& S_\text{PID, up} S_\text{up} + S_\text{PID, down} S_\text{down} \\
v_1 &=& \Psi_\text{up} \odot \Psi_\text{down} \odot \left(\frac{\partial q_\text{free, i}}{\partial h_\text{up, i}}\right)_{i = 1}^m + \text{diagonal}\left(J_I S_\text{PID, up}\right)\\
v_2 &=& q_\text{free} \odot \Psi_\text{down} \\
v_3 &=& \Psi_\text{up} \odot \Psi_\text{down} \odot \left(\frac{\partial q_\text{free, i}}{\partial h_\text{down, i}}\right)_{i = 1}^m + \text{diagonal}\left(J_I S_\text{PID, down}\right) \\
v_4 &=& q_\text{free} \odot \Psi_\text{up} \\
w_1 &=& h'(s) \\
w_2 &=& \psi'(s).
\end{align}
$$
This is highly structured, where $F$ denotes the local interaction between storages and flows, and the other term denotes the non-local influence via continuous control. For the number of Basins $n$ maxing out in the order of $10^5$ (which I think is comfortably the regime of Ribasim) and the number of continuous control nodes being very small w.r.t. $n$, Claude suggests a solve strategy via the Woodbury formula + sparse LU:
$$
\begin{align}
LU &=& I_n -\gamma MF \quad &\text{(sparse LU factorization)} \\
LUY &=& -\gamma M\frac{\partial q}{\partial c} \quad &\text{(apply factorization $\mu$ times in parallel)} \\
LUz_1 &=& c_s \quad &\text{(apply factorization once)} \\
\left(I_\mu + \frac{\partial c}{\partial s}Y\right) z_2 &=& \frac{\partial c}{\partial s}z_1 \quad &\text{($\mu$ sized dense solve)} \\
x_s &=& z_1 - Y z_2 \quad &\text{(explicit computation)}
\end{align}
$$
We can also exploit that the sparsity of $L$ and $U$ only has to be determined once for the whole simulation.
This is a bit out of date with #3113, but the general idea is still valid.
Let's investigate the structure of the linear system that emerges in #3113 for$x_s$ . The matrix there can be written as:
where
This is highly structured, where$F$ denotes the local interaction between storages and flows, and the other term denotes the non-local influence via continuous control. For the number of Basins $n$ maxing out in the order of $10^5$ (which I think is comfortably the regime of Ribasim) and the number of continuous control nodes being very small w.r.t. $n$ , Claude suggests a solve strategy via the Woodbury formula + sparse LU:
We can also exploit that the sparsity of$L$ and $U$ only has to be determined once for the whole simulation.