gwsource.diagnostics#

Budget residuals, closure and forward/backward duality checks, and the cycle report.

Exact identities of gross accounting (they follow from summing the unnormalised row equations
O_i F_i - sum_j q_ij F_j = d_i and I_i G_i - sum_j q_ji G_j = b_i over all cells; rows with
O_i = 0 or I_i = 0 satisfy them trivially). With r_i = I_i - O_i, Src_i and Snk_i the cell's
source and sink totals, d the destination sink flows and b the selected source flows:

    backward closure:  sum_i F_i Src_i - sum_i F_i r_i = sum_{k in D} |q_k|
    forward closure:   sum_i G_i Snk_i + sum_i G_i r_i = sum_{k in S} q_k
    duality:           sum_{k in S} F_{c(k)} q_k - sum_{k in D} G_{c(k)} |q_k| = sum_i F_i G_i r_i

So with the residual term the checks hold to round-off; without it they hold to within the
model's budget error (|sum F r| <= sum |r|, since 0 <= F <= 1).

Net accounting (the reference method [FB13]) uses intercell flows only: s_i = O^edges_i -
I^edges_i, T_i = I^edges_i + max(s_i, 0), f = 1 in the destination cells D (g = 1 in the source
cells S forward). There is no residual term, and the checks are

    net backward closure:  sum_i f_i max(s_i, 0)  = destination inflow
    net forward closure:   sum_i g_i max(-s_i, 0) = source inflow

where the net destination inflow is sum_D T_d - sum_{D->D} q - sum_{D->not D} q_ij f_j and the
net source inflow sum_S T_s - sum_{S->S} q - sum_{not S->S} q_ji g_j (gwsource.assemble). Both
sides use the solved fractions, so the net inflow depends on the solve, not on the budget alone,
and its meaning is pinned by tests/test_fb13_example.py and the hand-built tests. The error is
exactly the sum of the un-normalised row residuals over the cells outside the group,
sum_{i not in D} (T_i f_i - sum_j q_ij f_j) backward and sum_{i not in S} (T_i g_i - sum_j q_ji
g_j) forward, so the check tests the solve in sum.

Strong sinks (TraceConfig.strong_sinks). An absorbing cell a passes no water on, so its row
lacks the terms of its out-edges. Summing the rows as above leaves, in both accountings, the
strong-sink term A = sum_{a->j} q_aj (F_a - F_j) over the edges that leave absorbing cells: the
water they hand on physically but count with their own share F_a. The backward closure becomes

    sum_i F_i Src_i - sum_i F_i r_i - A = sum_{k in D} |q_k|      (gross)
    sum_i f_i max(s_i, 0) - A = destination inflow                 (net)

so the volume from no longer equals the destination's own inflow. With one group per sink term,
the fractions of every draining cell still sum to 1 (boundary control).

Cycle report: cycle_report() lists every strongly connected set of the flow graph with the flow
that goes round it (by cycle cancelling) and its share of the model throughflow; cycle_table()
and cycle_ids() turn it into rows and a per-cell array for export.
class gwsource.diagnostics.Check(lhs, rhs, residual_term, strong_sink_term=0.0)[source]#

Bases: object

Both sides of an identity and its error, with and without the residual term.
Parameters:
property error: float#
lhs - residual_term - strong_sink_term - rhs (zero to round-off).
property error_without_residual: float#
lhs - strong_sink_term - rhs, i.e. the error if the budget residual is ignored.
property relative_error: float#
|error| relative to |rhs| (or 1 if rhs is zero).
gwsource.diagnostics.cell_residual(step)[source]#
Return the per-cell budget residual r = I - O (positive: more water in than out).
Parameters:

step (FlowStep)

Return type:

ndarray

gwsource.diagnostics.budget_error(step)[source]#
Return the total absolute residual, total throughflow and their ratio.
Parameters:

step (FlowStep)

Return type:

dict[str, float]

gwsource.diagnostics.strong_sink_term(step, result, group)[source]#
Return A = sum q (F_a - F_j) over the edges leaving the group's absorbing cells (or 0).
Parameters:
Return type:

float

gwsource.diagnostics.backward_closure(step, result, group)[source]#
Return the closure check sum(volume from) - sum(F r) - A = destination inflow.

Net accounting uses intercell flows only, so its residual term is zero and the check is
sum(volume from) - A = destination inflow; A is zero without strong sinks (module docstring).
Parameters:
Return type:

Check

gwsource.diagnostics.forward_closure(step, result, group)[source]#
Return the closure check sum(volume to) + sum(G r) = source inflow.

Net accounting uses intercell flows only, so its residual term is zero and the check is
sum(volume to) = source inflow (module docstring).
Parameters:
Return type:

Check

gwsource.diagnostics.duality(step, back, dest_group, fwd, source_group)[source]#
Return the flow from a source group to a destination group, computed both ways (gross).
Parameters:
Return type:

Check

gwsource.diagnostics.cancel_cycles(edges, max_cancellations=100000)[source]#
Decompose the circulation of a flow graph by cycle cancelling.

Repeatedly find a directed cycle, subtract its smallest edge flow from every edge on it and
drop the emptied edge, until no cycle is left. Returns (circulation, circulating edge flow,
complete): circulation is the sum of the cancelled cycle flows; circulating edge flow is the
sum over cancelled cycles of cycle flow x cycle length, i.e. the throughflow the cells of the
set carry only because water goes round. The decomposition is not unique; it is a measure
of how much flow goes round, not a unique quantity. complete is False if max_cancellations
was reached first.
Parameters:
Return type:

tuple[float, float, bool]

gwsource.diagnostics.cycle_report(step, result=None)[source]#
Return the flow cycles of a step with their circulation, largest circulation first.

For each strongly connected set: the cells, the flow on its internal edges, the circulation
found by cycle cancelling, the circulating throughflow (circulation x cycle length) and its
share of the total model throughflow (sum of cell outflows). Sets reported as closed loops
by a trace result are flagged closed.
Parameters:
Return type:

list[Cycle]

gwsource.diagnostics.cycle_table(cycles)[source]#
Return one row (dict) per cycle: id, size, flows, share and the 0-based cells.
Parameters:

cycles (list[Cycle])

Return type:

list[dict]

gwsource.diagnostics.cycle_ids(n_cells, cycles)[source]#
Return, per cell, the index of the cycle it belongs to, or -1.
Parameters:
Return type:

ndarray