Solver divergence with large quadratic resistances.

[natalia-rubio]

Description

When a 0D system features large stenosis (quadratic) resistances, especially if they vary in sign. The solver often diverges, that is, the linear Newton solver in the time-stepping scheme is not able to sufficiently reduce the residual. This is somewhat expected, because the introduction of large quadratic resistances of varying sign means the relationship between the solution variables and the residual is not monotonic or convex. There is technically no guarantee that a solution exists.

Reproduction

bifurcations_EL_NN_JunctionOnly.json

Sample input files with large quadratic resistances that diverges.

Expected behavior

Ideally, the solver should converge if there is a solution, and print an error if not.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[zasexton]

working on this in branch issue/237-large-quadratic-resistance-divergence.

the current branch adds robust Newton controls for stiff cases, including damped Newton line search, residual/variable scaling, adaptive timestep retries. Initial testing shows these changes improve robustness on the clamped large bifurcation case (when R, L, and S >= 0), and we are continuing to evaluate performance and convergence behavior to determine in which parameter regimes solution envelopes exist.

[zasexton]

Okay after reviewing the input parameter conditions with sympy symbolic analysis over the last few days I have the following general thoughts for blood vessel parameters:

The steady problem that is initially posed on a undirected network with the incidence/constraint matrix $B$: find edge flows $q$ and nodal pressures $p$ satisfying

$$Bq = b, \qquad B^\top p = \phi(q), \qquad \phi_e(q_e) = R_e q_e + S_e |q_e| q_e,$$

with incremental resistance $\phi_e'(q_e) = R_e + 2 S_e |q_e|$. This system is exactly the KKT (stationarity) condition of the constrained minimization

$$\min_q ; \sum_e \psi_e(q_e) \quad \text{s.t.} \quad Bq = b, \qquad \psi_e(q_e) = \tfrac{1}{2} R_e q_e^2 + \tfrac{1}{3} S_e |q_e|^3,$$

with the pressures as Lagrange multipliers of the conservation constraints. Two structural conditions are required independent of the edge coefficients:

1. net inflows sum to zero (conservation of mass) on every connected component lacking a fixed-pressure reference
2. one grounded node per connected component to stipulate the pressure gauge

Solvability is then governed by convexity and coercivity of the constrained energy program. As background, passivity means that the network components do not create any energy internally (i.e. increases in energy should come from boundary/input power). Passivity can be determined from the instantaneous power of vessels within a network. Since the pressure drop in the steady case can be given as $\Delta p_e = R_e + S_e |q_e|q_e$, the instantaneous power dissipated by a given edge is $q_e \Delta p = R_e q_e^2 + S_e |q_e|^3 = |q_e|^2 (R_e + S_e |q_e|)$. Therefore a passive network will satisfy $|q_e| <= R_e/(-S_e)$ and means that the network stored energy can only increase due to boundary nodes of the system. This is a weaker condition than uniqueness but it is important because it can lead to nonphysical growth that may cause instability for the transient solve.

With $R_e \ge 0$ and $S_e \ge 0$ on every edge, each $\psi_e$ is convex and each $\phi_e$ is monotone. If $S_e < 0$ is admitted, only envelope-local guarantees exist. Within a certified bound $|q_e| \le Q_e$:

• Uniqueness and Newton robustness require a local type of incremental monotonicity $R_e + 2 S_e Q_e > 0$, i.e. $Q_e < R_e / (-2 S_e)$ — a factor of two stricter than the passivity condition $R_e + S_e Q_e \ge 0$.

• The network-level certificate is $Z^\top \mathrm{diag}(d_e), Z \succ 0$ with worst-case diagonals $d_e = R_e + 2 S_e Q_e$ for $S_e < 0$ and $d_e = R_e$ otherwise, where $Z$ spans the feasible null space.

• Because $\psi_e \to -\infty$ when $S_e < 0$, no global minimizer exists and solutions may lie outside the envelope; $Q_e$ must be certified a priori.

So the hierarchy is:

$R + 2S|q| > 0$ => monotone, unique/Newton-friendly
$R + S|q| >= 0$ => passive/dissipative
$R + S|q| < 0$ => active/energy-generating

There is an intermediate region for $S < 0$:

$R / (-2S) < |q| <= R / (-S)$

So while we want to remain in the unique/Newton-friendly regime, we always want to stay physically plausible. Being in the physically plausible regime does not necessarily guarantee that Newton will succeed.

Edges with $R_e = 0$, $S_e > 0$ have $\phi_e'(0) = 0$ — which suggests a degeneracy that stalls Newton near zero flow — so input files should probably enforce $R_e \ge \varepsilon > 0$ to maintain good Newton performance (also it would seem physically implausible for a section of vessels to have absolutely no nominal resistance).

When the solver proceeds to the transient portion of the simulation. The dynamics enter through the capacitance and inductance:

$$H(p, q) = \tfrac{1}{2}, p^\top C, p + \tfrac{1}{2}, q^\top L, q,$$

and passive storage requires $C \ge 0$ and $L \ge 0$, preferably strictly positive on all independent dynamic states. With $C \succ 0$, $L \succ 0$, $R_e \ge 0$, $S_e \ge 0$, the system is a dissipative ODE with $dH/dt$ bounded by boundary power and implicit stepping is strictly monotone because the step diagonal $L_e/\Delta t + R_e + 2 S_e |q_e|$ is positive. Zeros in $C$ or $L$ degrade the ODE to a DAE whose index depends on topology: compliance-free nodes make their pressures algebraic, and lossless loop configurations may necessitate that initial conditions cannot be arbitrary and must also satisfy algebraic constraints for state variables $p_e(0)$, $q_e(0)$ and potentially derivatives $p'_e(0)$, $q'_e(0)$. Negative incremental resistance can be tolerated transiently only through the timestep-dependent bound

$$\Delta t < \frac{L_e}{-R_e - 2 S_e Q_e} \quad \text{whenever } R_e + 2 S_e Q_e < 0,$$

but this is numerical regularization, not passivity: such operating points are generically unstable and adaptive step-size growth can silently cross the regularity threshold mid-run.