gwsource.solve#

Solve (I - P) X = B for many right-hand sides; detect, report and (net mode) remove cycles.

Solvers:
- topological: if the dependency graph of the matrix is acyclic, permute it to lower-triangular
  form (dependencies first) and solve by substitution (scipy's compiled spsolve_triangular).
  This is the reference method's algorithm [FB13] in matrix form.
- block: permute the matrix to block lower-triangular form, one diagonal block per strongly
  connected component in topological order, and factorise it with SuperLU in that natural
  order (no fill-reducing column permutation, no pivoting). Fill-in is then confined to the
  rows and columns of the cyclic blocks, so a 1M-cell graph with thousands of small cycles
  factorises in well under a second. Works for any graph; the default for cyclic graphs.
- lu: sparse LU with SuperLU's default COLAMD ordering on the original matrix, factorised once
  and reused for every column. Kept as a reference; its fill-in makes it impractical beyond
  about 1e5 cells of a 3D grid.
- auto: topological if the graph is acyclic, otherwise block.

factorize() returns the factorisation (a Factor) so that further right-hand sides that depend on
the first solution (the travel-time systems of gwsource.traveltime) reuse it.

Pivoting is not needed: I - P is a (weakly) row-diagonally dominant M-matrix. It is singular only
on a closed recirculation loop, a strongly connected set with no exit: no dependency leaves it
and no row in it has an exit share 1 - sum_j P_ij (a sink, or in forward mode a source) above a
tolerance. isolate_closed_loops() finds these before the solve, turns their rows into identity
rows with a zero right-hand side (fraction 0: water circulating in them never reaches, or never
came from, any group) and reports them.

Cycles are the strongly connected components with more than one cell of the positive-flow
digraph (scipy.sparse.csgraph.connected_components with connection="strong"). scipy labels the
components in reverse topological order of the condensed graph; the order is verified, with
Kahn's algorithm as a fallback.
class gwsource.solve.Cycle(cells, internal_flow, min_edge_flow, circulation=nan, circulating_throughflow=nan, throughflow_share=nan, closed=False)[source]#

Bases: object

One strongly connected set of cells in the flow graph (a flow cycle or network of cycles).

internal_flow is the sum of the flows on the edges inside the set; min_edge_flow the
smallest of them. The circulation fields are filled by diagnostics.cycle_report().
Parameters:
property n_cells: int#
Number of cells in the cycle set.
gwsource.solve.strong_components(n, src, dst)[source]#
Return (number of components, labels) of the strongly connected components of a digraph.
Parameters:
Return type:

tuple[int, ndarray]

gwsource.solve.find_cycles(n, edges)[source]#
Return every strongly connected set of more than one cell, with its internal edge flows.
Parameters:
Return type:

list[Cycle]

gwsource.solve.remove_cycle_edges(n, edges)[source]#
Return a keep-mask deleting the smallest edge of each cycle until the graph is acyclic.
Parameters:
Return type:

ndarray

gwsource.solve.component_order(n, before, after, labels)[source]#
Return a cell order with the components in topological order (cells of one together).
Parameters:
Return type:

ndarray

gwsource.solve.topological_order(n, before, after)[source]#
Return an order in which every 'before' cell precedes its 'after' cell, or None if cyclic.
Parameters:
Return type:

ndarray | None

gwsource.solve.closed_components(matrix, labels, tol=1e-12)[source]#
Return a cell mask of the strongly connected sets (> 1 cell) of I - P with no exit.
Parameters:
Return type:

ndarray

class gwsource.solve.DependencyGraph(before, after, n_components, labels)[source]#

Bases: object

Dependency edges of a system matrix and their strongly connected components.
Parameters:
gwsource.solve.dependency_graph(matrix)[source]#
Return the dependency edges of matrix and their strongly connected components.
Parameters:

matrix (csr_matrix)

Return type:

DependencyGraph

gwsource.solve.isolate_closed_loops(matrix, rhs, tol=1e-12, graph=None)[source]#
Return (matrix, rhs, loops) with closed recirculation loops set to identity rows, rhs 0.
Parameters:
Return type:

tuple[csr_matrix, ndarray, list[ndarray]]

gwsource.solve.identity_rows(matrix, closed)[source]#
Return matrix with the off-diagonal entries of the closed rows removed.
Parameters:
Return type:

csr_matrix

class gwsource.solve.Factor(kind, order=None, triangular=None, lu=None, zero_rows=None)[source]#

Bases: object

A factorised system matrix, reusable for further right-hand sides (e.g. travel times).

kind is the solver used; order the cell order of the permuted matrix (None for "lu");
triangular the permuted lower-triangular matrix ("topological"); lu the SuperLU object
("block" and "lu"). zero_rows marks closed-loop cells whose rows were made identity rows:
their right-hand side is forced to 0 by solve(), as isolate_closed_loops() did for the first
solve. Not picklable (SuperLU objects are not); keep it within one process.
Parameters:
solve(rhs)[source]#
Return X with matrix @ X = rhs (rows of closed loops forced to 0).
Parameters:

rhs (ndarray)

Return type:

ndarray

gwsource.solve.factor_topological(matrix, order)[source]#
Permute the matrix to lower-triangular form for substitution.
Parameters:
Return type:

Factor

gwsource.solve.factor_lu(matrix)[source]#
Factorise with SuperLU and its default COLAMD ordering.
Parameters:

matrix (csr_matrix)

Return type:

Factor

gwsource.solve.factor_block(matrix, order)[source]#
Factorise with SuperLU in block lower-triangular (component) order, without reordering.
Parameters:
Return type:

Factor

gwsource.solve.factorize(matrix, solver='auto', graph=None)[source]#
Return a reusable factorisation of matrix with the chosen solver (see module docstring).

graph, if given, must be dependency_graph(matrix) (it saves recomputing the components).
Parameters:
Return type:

Factor

gwsource.solve.solve_topological(matrix, rhs, order)[source]#
Solve by substitution after permuting the matrix to lower-triangular form.
Parameters:
Return type:

ndarray

gwsource.solve.solve_lu(matrix, rhs)[source]#
Solve with one sparse LU factorisation reused for every right-hand-side column.
Parameters:
Return type:

ndarray

gwsource.solve.solve_block(matrix, rhs, order)[source]#
Solve by SuperLU in block lower-triangular (component) order, without reordering.
Parameters:
Return type:

ndarray

gwsource.solve.solve_system(matrix, rhs, solver='auto', graph=None)[source]#
Solve matrix @ X = rhs; return (X, name of the solver used).

graph, if given, must be dependency_graph(matrix) (it saves recomputing the components).
Parameters:
Return type:

tuple[ndarray, str]