gwsource.experimental.face_routing#

Face-to-face routing in rectilinear DIS cells (Chowdhury 2016), an alternative to full mixing.

Full mixing (gwsource.trace) assumes that water entering a cell leaves through every outflow face
in proportion to that face's flow, whatever face it came in through. Here each inflow face's
water is partitioned among the outflow faces with Pollock's (1988) semi-analytical velocity
field, as proposed by Chowdhury (2016), and the trace runs on a graph whose nodes are the cell
faces (connections) rather than the cells.

Velocity field (Pollock 1988)
-----------------------------
In a cell of size dx x dy each velocity component varies linearly along its own axis only:
v_x(x) between the face velocities v_x1 and v_x2 (face flow / face area), v_y(y) likewise. A
particle's motions along x and along y are independent, and the time to reach a face along one
axis depends only on the start position s along that axis:
    t = ln(v_face / v(s)) / A,   A = (v2 - v1) / L      (t = distance / v when A = 0),
infinite if the particle heads for an inflow or no-flow face (a stagnation point lies between).
A particle leaves through the face it reaches first.

Partition of a lateral inflow face. Water enters face x1 (say) spread uniformly over the face
(Pollock's assumption). Every particle needs the same time T to cross the cell along x (infinite
if x2 is not an outflow face); along y its exit time depends only on its entry position y0. The
positions whose particles reach face y2 before T form an interval of closed-form length
    m_y2 = min(v_y2 (1 - exp(-A_y T)) / (A_y dy), share of positions moving towards y2)
(m_y1 is the mirror image); the remainder 1 - m_y1 - m_y2 leaves through x2. This reproduces
Chowdhury's (2016) case formulas for cells with two inflow and two outflow faces (stagnation-point
and corner cases) and also covers the cells with one inflow or one outflow face.

Areal sources: Pollock's treatment of distributed recharge
----------------------------------------------------------
A source term can be put on a cell face with a MODPATH IFACE code (FaceRoutingConfig.iface; by
default recharge enters the top face, IFACE 6). A routed cell may have flow through one vertical
face only, and it must be an inflow (recharge on top of a single-layer model, or leakage from the
layer above): the vertical velocity then falls linearly to 0 at the opposite, closed face, so
water entering through a lateral face never reaches a vertical face and the lateral rows above
are unchanged. The vertical inflow takes what is left of each outflow face by mass balance,
    W[top, g] proportional to max(Q_g - sum_f Q_f W[f, g], 0),
which is exact for Pollock's field (its face flows balance).

Which cells are routed
----------------------
A cell is routed if every external term in it with non-zero flow is assigned to a face, no face
carries both inflow and outflow, at most one vertical face has flow and it is an inflow, and the
cell has both inflow and outflow faces. Every other cell is fully mixed exactly as in
gwsource.trace: cells with internal (IFACE 0) sources or sinks such as wells, river or drain
leakage and constant heads; cells with flow through both vertical faces (most cells of a
multi-layer model: this prototype is 2D); strong sinks and sources (where routing and mixing
coincide). Mixing those cells is a choice, not a derivation: water passing a weak sink is not
split by the face it entered through, so weak-sink cells keep the full-mixing error.

Graph
-----
Unknowns: one value per cell (used by mixed cells; routed cells get identity rows) and one per
positive-flow connection, the "face node" e (the water crossing that face). For a face node e
entering cell j through face f:
    routed j:  F_e = sum_g W_j[f, g] (sum_{b on g} q_b / Q_g F_b + sum_{sinks s on g} q_s / Q_g d_s)
    mixed j:   F_e = F_j,
and a mixed cell is a cell node of gross accounting: F_j = sum_e q_e / O_j F_e + sum_s q_s / O_j d_s
(d_s = 1 for a sink term of the destination group). Solved with gwsource.solve (one
factorisation, one column per group). A source term on face f of a routed cell takes
F_k = sum_g W_j[f, g] (...), a source term in a mixed cell F_k = F_j.

Per-cell outputs (comparable with gwsource.trace): volume_through_j = the sum of q F over the
cell's outflow connections plus its destination sink flows; fraction_through = that / O_j (the
outflow-weighted mean over the cell's outflow faces; F_j itself in mixed cells);
volume_from = the sum over the cell's source terms of q_k F_k.

Closure. With In_m the water a node receives through the plans, r_m = In_m - O_m (mixed cells: the
cell residual; face nodes and sinks of a routed cell: their share of that cell's residual, since
rows of W sum to 1 but its columns match the outflow faces only to the flow model's precision):
    sum_k q_k F_k - sum_m r_m F_m = destination inflow
exactly (to round-off), as in gwsource.diagnostics.backward_closure.

Limitations: DIS only; 2D (at most one vertical inflow face per routed cell); steady flow (one
FlowStep); backward mode only; no travel time. Water is still mixed across each face (where a
streamline crosses a face is not kept), so flow oblique to the grid keeps some lateral spreading,
much less than full mixing.

IFACE codes (MODPATH): 0 internal, 1 column-1 side (x1), 2 column+1 side (x2), 3 row+1 side (y1),
4 row-1 side (y2), 5 bottom (layer+1 side), 6 top (layer-1 side). Arrays here use face index
IFACE - 1 (0..5), and -1 for internal.
class gwsource.experimental.face_routing.FaceRoutingConfig(iface=<factory>, route=True, solver='auto', closed_loop_tol=1e-12)[source]#

Bases: object

Options of route_backward().

iface: {package name or budget text: MODPATH IFACE code} for terms that enter or leave
    through a cell face (default: recharge through the top face); all other terms are
    internal and make their cell fully mixed.
route: False mixes every cell (reproduces gwsource.trace.backward on the face graph).
solver: as TraceConfig.solver.
closed_loop_tol: as TraceConfig.closed_loop_tol.
Parameters:
class gwsource.experimental.face_routing.CellFaces(outflow, inflow)[source]#

Bases: object

Per-cell (n, 6) gross outflow and inflow through each face (connections + face terms).
Parameters:
property net: ndarray#
Net outflow per face (positive out of the cell).
class gwsource.experimental.face_routing.FaceGraph(edges, out_face, in_face, routed, plans, transition, sink_share, source_map, source_sink, residual, sink_residual)[source]#

Bases: object

The face-routing system: nodes are cells (0..n-1) then connections (n + edge index).

transition: P (nodes x nodes); sink_share: nodes x terms; source_map / source_sink: terms x
nodes / terms x terms (F of a source term = source_map X + source_sink D); residual: per node;
sink_residual: per term (sinks of routed cells).
Parameters:
class gwsource.experimental.face_routing.FaceRoutingResult(selections, graph, fraction_through, volume_through, volume_from, volume_from_by_term, destination_inflow, face_fraction, term_fraction, residual_term, solver_used='', closed_loops=<factory>)[source]#

Bases: object

Backward face-routing outputs per destination group (per-cell arrays as in BackwardResult).

face_fraction[g]: F of each face node (graph.edges order); term_fraction[g]: F of each
external term (meaningful for source terms); residual_term[g]: the closure residual term.
Parameters:
property routed: ndarray#
Cell mask of the routed (not fully mixed) cells.
gwsource.experimental.face_routing.read_dis_spacing(grb_file)[source]#
Return (delr, delc) from a DIS binary grid file.
Parameters:

grb_file (str)

Return type:

tuple[ndarray, ndarray]

gwsource.experimental.face_routing.check_dis(step, delr, delc)[source]#
Assert that the step is a DIS grid matching delr/delc; return (nlay, nrow, ncol).
Parameters:
Return type:

tuple[int, int, int]

gwsource.experimental.face_routing.neighbour_face(src, dst, shape)[source]#
Return the face index (IFACE - 1) of cell src through which it connects to dst.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.term_faces(step, iface)[source]#
Return the face index of every external term (-1 for internal terms).
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.cell_faces(n, edges, out_face, step, term_face)[source]#
Return per-cell outflow and inflow through each face (connections + face-assigned terms).
Parameters:
Return type:

CellFaces

gwsource.experimental.face_routing.routable_cells(step, faces, term_face)[source]#
Return the cells that can be routed (see module docstring); all others are mixed.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.crossing_time(v_in, v_out, length)[source]#
Return the time to cross a cell from an inflow face (speed v_in > 0) to the opposite face.

v_out is the velocity at the opposite face, positive onwards; inf if it is not an outflow.
Parameters:
Return type:

float

gwsource.experimental.face_routing.exit_measures(v1, v2, length, t)[source]#
Return the shares of uniformly spread start positions along an axis leaving through face
1 and face 2 within time t (v1, v2: face velocities in the +axis direction).
Parameters:
Return type:

tuple[float, float]

gwsource.experimental.face_routing.lateral_row(face, net, dx, dy)[source]#
Return the shares of a lateral inflow face's water leaving through each face.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.vertical_row(plan, net)[source]#
Return the vertical inflow's shares: what the lateral inflows leave of each outflow face.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.cell_plan(net, dx, dy)[source]#
Return W[inflow face, outflow face]: share of each inflow face's water per outflow face.

net: (6,) net outflow per face (IFACE order); rows of non-inflow faces are zero.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.build_graph(step, delr, delc, config)[source]#
Return the face-routing graph of a DIS step (see module docstring).
Parameters:
Return type:

FaceGraph

gwsource.experimental.face_routing.destination_matrix(step, groups)[source]#
Return D (terms x groups): 1 where a selected term is a sink.
Parameters:
Return type:

ndarray

gwsource.experimental.face_routing.route_backward(step, groups, delr, delc, config=None)[source]#
Trace backward from each destination group with face-to-face routing in DIS cells.
Parameters:
Return type:

FaceRoutingResult

gwsource.experimental.face_routing.closure(result, group)[source]#
Return the closure check sum(volume from) - residual term = destination inflow.
Parameters:
Return type:

Check