Release 1.4.0: PDN mode, the config-file workflow, and the dialog editor

Multiple Thevenin supplies and prescribed-current loads on one net,
solved in absolute volts with the Tellegen power balance verified per
run; a source-sink pair table (effective copper resistance per
supply x load pair plus an exactly-summing proportional-sharing loss
attribution), in summary.txt and as its own figure. Bonded terminals
short a package's contacts into one lug so the per-pin split becomes
a solve outcome. Geometry dumps carry the terminal set (schema v8).

The dialog gained a Classic/PDN mode selector and a full PDN editor:
per-role supply/load tables built from the marker rectangles (or a
config's terminal set, which never pins mode or net), with Component
hints, per-terminal Layer scopes, Active checkboxes, comments, a
per-net row filter, resizable tables and a scrolling, screen-sized
dialog. Numbers accept SI suffixes (50m, 4.7k) everywhere.

fill_res_config.json fully specifies a run (classic or PDN) with
validation, comments, named side-by-side configs (the one called
default auto-loads), Load/Save buttons with an editable file name,
and saves that never drop anything drawn on the board.

347 tests, green on Python 3.13 and on the 3.9 macOS wheel stack.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

fill_resistance/solver.py

@@ -59,12 +59,19 @@ class SolveInfo:
     residual: float | None = None
 
 
+# via_index tag for PDN supply-attachment edges (virtual Thevenin node
+# to contact cell): excluded from the in-plane fields (== -1) AND from
+# the via power/reports (>= 0); their dissipation is P_supply_internal
+PDN_EDGE = -2
+
+
 @dataclass
 class Edges:
     a: np.ndarray                 # int64 flat cell ids
     b: np.ndarray
     w: np.ndarray                 # conductance [S]
-    via_index: np.ndarray         # int32; -1 = in-plane edge
+    via_index: np.ndarray         # int32; -1 = in-plane edge,
+                                  # PDN_EDGE = supply attachment
     dead_barrels: int = 0         # barrels spanning >=2 layers that found
                                   # fill copper on fewer than 2 of them
 
@@ -79,6 +86,59 @@ class ViaReport:
     power_w: float                # total barrel dissipation @ I_test
 
 
+@dataclass
+class SupplyReport:
+    """One PDN supply after the solve. The delivered current is an
+    OUTCOME (Thevenin split), not an input."""
+    label: str
+    v_oc: float                   # open-circuit volts used in the solve
+    r_out_ohm: float
+    i_a: float                    # delivered current [A]
+    v_contact: float              # mean volts over the contact cells
+    p_internal_w: float           # dissipated inside r_out
+    part_currents: list = field(default_factory=list)  # [(label, amps)]
+    v_eff: float = 0.0            # current-weighted contact volts: the
+                                  # potential the delivered power sees
+                                  # (= v_contact for ideal and bonded
+                                  # contacts); makes the pair-loss
+                                  # allocation sum EXACTLY to the
+                                  # copper dissipation
+    component: str = ""           # display: terminal's owner hint
+    comment: str = ""             # display: terminal's free-text note
+
+
+@dataclass
+class LoadReport:
+    """One PDN load after the solve. The draw is prescribed; the contact
+    voltage is the outcome of interest."""
+    label: str
+    i_a: float                    # prescribed draw [A]
+    v_mean: float                 # mean volts over the contact cells
+    v_min: float                  # worst-case contact cell
+    p_w: float                    # i_a * v_mean (exact: injection and
+                                  # averaging weights coincide)
+    part_currents: list = field(default_factory=list)  # [(label, amps)]
+    component: str = ""           # display: terminal's owner hint
+    comment: str = ""             # display: terminal's free-text note
+
+
+@dataclass
+class PairReport:
+    """One (supply, load) pair: the effective COPPER resistance between
+    the two contacts (source internals excluded; injection patterns as
+    in the solve - uniform per cell, or the bonded lug) and the copper
+    loss attributed to the pair by PROPORTIONAL SHARING
+    (f_ij = I_i * I_j / I_component, P_ij = f_ij * (v_eff_i - v_mean_j)).
+    The attribution is a convention, not unique physics - but it sums
+    exactly to the total copper dissipation, and R is an operating-
+    point-independent property of the board."""
+    supply: str
+    load: str
+    r_ohm: float | None           # None: no common copper path
+    i_share_a: float              # attributed current [A]
+    p_w: float                    # attributed copper loss [W]
+
+
 @dataclass
 class Result:
     R_ohm: float
@@ -106,6 +166,18 @@ class Result:
     skin_depth_um: float | None = None
     rs_ratios: list[float] = field(default_factory=list)  # R_AC/R_DC per layer
     timings: dict = field(default_factory=dict)
+    # --- PDN mode (mode == "pdn"; classic solves leave these empty) ---
+    # R_ohm is NaN there (no single two-terminal R); i_test carries the
+    # summed load draw so %-of-total displays keep working; fields V/
+    # Jmag/Parea are in ABSOLUTE volts / real operating current
+    mode: str = "classic"                     # "classic" | "pdn"
+    supplies: list = field(default_factory=list)   # [SupplyReport]
+    loads: list = field(default_factory=list)      # [LoadReport]
+    P_loads: float = 0.0                      # sum of load powers [W]
+    P_supply_internal: float = 0.0            # sum of r_out dissipation
+    v_nominal: float | None = None            # default supply v_oc used
+    pairs: list = field(default_factory=list)  # [PairReport], every
+                                               # supply x load
 
 
 def _shifts2d():
@@ -288,20 +360,116 @@ def connected_restrict(stack: RasterStack, e1: np.ndarray, e2: np.ndarray,
     return changed, len(common)
 
 
-def _assemble(state: np.ndarray, edges: Edges, rhs_extra: np.ndarray | None):
+def _pdn_keep_components(terminals: list, per_term: list) -> set:
+    """The PDN component keep rule + its diagnostics, on the label sets
+    each terminal's contact touches (shared by the cell graph here and
+    the adaptive leaf graph). Keep components holding >= 1 supply AND
+    (>= 1 load OR >= 2 supplies); errors for unreachable / sheet-
+    spanning loads, notes for sheet-spanning supplies."""
+    n_sup: dict = {}
+    has_load: set = set()
+    for t, labs in zip(terminals, per_term):
+        if t.role == "supply":
+            for l in labs:
+                n_sup[l] = n_sup.get(l, 0) + 1
+        else:
+            has_load |= labs
+    kept = {l for l, c in n_sup.items() if l in has_load or c >= 2}
+    if not kept:
+        raise ConnectivityError(
+            "No copper component connects a supply to a load (not even "
+            "through vias). Check the layer selection, the terminal "
+            "definitions and that the fills are up to date."
+        )
+    for t, labs in zip(terminals, per_term):
+        if t.role != "load":
+            continue
+        kl = labs & kept
+        if not kl:
+            raise ConnectivityError(
+                f"Load '{t.label}' sits on copper that is not connected "
+                f"to any supply (not even through vias)."
+            )
+        if len(kl) > 1:
+            if t.bonded:
+                # the external bond IS the connection: the split
+                # between the sheets is well-defined through the lug
+                print(f"note: bonded load '{t.label}' spans {len(kl)} "
+                      f"disconnected copper sheets; the split between "
+                      f"them is set by its external bond")
+                continue
+            raise ConnectivityError(
+                f"Load '{t.label}' spans {len(kl)} disconnected copper "
+                f"sheets - the current split between them is undefined "
+                f"with per-cell injection. Include the layers/vias that "
+                f"join them, mark the load as bonded, or split it into "
+                f"one terminal per sheet."
+            )
+    for t, labs in zip(terminals, per_term):
+        if t.role == "supply" and len(labs & kept) > 1:
+            print(f"note: supply '{t.label}' feeds {len(labs & kept)} "
+                  f"disconnected copper sheets; the split between them "
+                  f"is set by its output resistance (Thevenin)")
+    return kept
+
+
+def connected_restrict_multi(stack: RasterStack, term_masks: list,
+                             terminals: list, edges: Edges
+                             ) -> tuple[bool, int]:
+    """PDN connectivity restriction: keep copper components holding
+    >= 1 supply AND (>= 1 load OR >= 2 supplies) - the second clause
+    keeps circulating-current paths between paralleled supplies with
+    unequal v_oc. A load on copper reachable from no supply is an
+    error; so is a load spanning several kept components (its uniform
+    injection cannot decide the split between disconnected sheets). A
+    load merely LOSING cells to dropped copper is fine: those cells
+    could not carry current anyway, the draw renormalizes over the
+    rest. Mutates stack.masks and the term_masks. Returns (changed,
+    n_kept_components)."""
+    n = stack.masks.size
+    graph = sparse.coo_matrix(
+        (np.ones(len(edges.a)), (edges.a, edges.b)), shape=(n, n))
+    _, labels = csgraph.connected_components(graph, directed=False)
+    labels3 = labels.reshape(stack.masks.shape)
+
+    per_term = [set(np.unique(labels3[m]).tolist()) if m.any() else set()
+                for m in term_masks]
+    kept = _pdn_keep_components(terminals, per_term)
+    keep = np.isin(labels3, sorted(kept)) & stack.masks
+    changed = bool((stack.masks & ~keep).any())
+    stack.masks &= keep
+    for t, m in zip(terminals, term_masks):
+        had = bool(m.any())
+        m &= keep
+        if t.role == "supply" and had and not m.any():
+            print(f"warning: supply '{t.label}' only touches copper not "
+                  f"connected to any load - it delivers 0 A")
+    return changed, len(kept)
+
+
+def _assemble(state: np.ndarray, edges: Edges, rhs_extra: np.ndarray | None,
+              dirichlet_v: np.ndarray | None = None):
     """Weighted-Laplacian assembly with Dirichlet elimination.
     state: 0 off, 1 free, 2 Dirichlet@1V, 3 Dirichlet@0V.
-    rhs_extra: per-flat-cell current injection [A] added for free cells."""
+    rhs_extra: per-flat-cell current injection [A] added for free cells.
+    dirichlet_v (PDN mode): per-node Dirichlet volts - any state >= 2
+    is then held at dirichlet_v[node] instead of the fixed 1 V / 0 V
+    pair, and the direct-connection short check is skipped (edges
+    between Dirichlet nodes simply conduct; their currents come out of
+    the post-solve edge fluxes). None keeps the classic behavior
+    bit-for-bit."""
     n = state.size
     sa, sb = state[edges.a], state[edges.b]
-    short = ((sa == 2) & (sb == 3)) | ((sa == 3) & (sb == 2))
-    if short.any():
-        n_via = int((edges.via_index[short] >= 0).sum())
-        raise ElectrodeError(
-            f"The terminals are directly connected by {int(short.sum())} "
-            f"conductance(s) ({n_via} via barrel(s)) without any free copper "
-            f"in between - move the contacts apart."
-        )
+    if dirichlet_v is None:
+        short = ((sa == 2) & (sb == 3)) | ((sa == 3) & (sb == 2))
+        if short.any():
+            n_via = int((edges.via_index[short] >= 0).sum())
+            raise ElectrodeError(
+                f"The terminals are directly connected by "
+                f"{int(short.sum())} conductance(s) ({n_via} via "
+                f"barrel(s)) without any free copper in between - move "
+                f"the contacts apart."
+            )
 
     free = state == 1
     n_free = int(free.sum())
@@ -318,10 +486,18 @@ def _assemble(state: np.ndarray, edges: Edges, rhs_extra: np.ndarray | None):
     fa, fb = sa == 1, sb == 1
     np.add.at(diag, idx[edges.a[fa]], edges.w[fa])
     np.add.at(diag, idx[edges.b[fb]], edges.w[fb])
-    r1a = fa & (sb == 2)
-    r1b = fb & (sa == 2)
-    np.add.at(rhs, idx[edges.a[r1a]], edges.w[r1a])
-    np.add.at(rhs, idx[edges.b[r1b]], edges.w[r1b])
+    if dirichlet_v is None:
+        r1a = fa & (sb == 2)
+        r1b = fb & (sa == 2)
+        np.add.at(rhs, idx[edges.a[r1a]], edges.w[r1a])
+        np.add.at(rhs, idx[edges.b[r1b]], edges.w[r1b])
+    else:
+        r1a = fa & (sb >= 2)
+        r1b = fb & (sa >= 2)
+        np.add.at(rhs, idx[edges.a[r1a]],
+                  edges.w[r1a] * dirichlet_v[edges.b[r1a]])
+        np.add.at(rhs, idx[edges.b[r1b]],
+                  edges.w[r1b] * dirichlet_v[edges.a[r1b]])
     if rhs_extra is not None:
         rhs += rhs_extra[free]
 
@@ -605,6 +781,71 @@ def _part_currents(parts, Ie, edges, e_flat, scale,
     return out
 
 
+def _postprocess_fields(problem: Problem, stack: RasterStack, edges: Edges,
+                        Vflat: np.ndarray, s: float, sigmas: list[float],
+                        sigma_buildup: float):
+    """Edge powers, per-layer dissipation, via reports and the V/J/P
+    display fields on the uniform grid - shared verbatim by the classic
+    and PDN solves. PDN appends virtual supply nodes after the grid
+    ids: everything here slices Vflat back to the grid (a no-op view
+    for classic) and selects in-plane edges by via_index == -1 / via
+    barrels by >= 0, so PDN_EDGE attachment edges stay out of the
+    copper fields and the via reports. Returns (Pe, Ie, P_layers,
+    P_vias, Parea, via_reports, V3, J3): Pe is s^2-scaled, Ie is at the
+    drive of Vflat (unit drive classic, absolute volts PDN); via report
+    currents are s-scaled here."""
+    L, ny, nx = stack.masks.shape
+    h_m = stack.h_nm * 1e-9
+    n_grid = stack.masks.size
+
+    # per-edge power @ I_test; distribute in-plane power to endpoint cells
+    Pe = edges.w * ((Vflat[edges.a] - Vflat[edges.b]) * s) ** 2
+    inplane = edges.via_index == -1
+    Pflat = np.zeros(n_grid)
+    np.add.at(Pflat, edges.a[inplane], 0.5 * Pe[inplane])
+    np.add.at(Pflat, edges.b[inplane], 0.5 * Pe[inplane])
+    Parea = Pflat.reshape(L, ny, nx) / (h_m * h_m)
+    Parea[~stack.masks] = np.nan
+    plane = ny * nx
+    P_layers = [float(Pflat[li * plane:(li + 1) * plane].sum())
+                for li in range(L)]
+    P_vias = float(Pe[edges.via_index >= 0].sum())
+
+    # via reports: max segment current + total power per via
+    Ie = edges.w * (Vflat[edges.a] - Vflat[edges.b])
+    via_reports = []
+    if problem.vias:
+        vidx = edges.via_index
+        for vi in np.unique(vidx[vidx >= 0]):
+            sel = vidx == vi
+            via = problem.vias[vi]
+            via_reports.append(ViaReport(
+                x_mm=via.x * 1e-6, y_mm=via.y * 1e-6, kind=via.kind,
+                drill_mm=via.drill_nm * 1e-6,
+                current_a=float(np.abs(Ie[sel]).max()) * s,
+                power_w=float(Pe[sel].sum()),
+            ))
+        via_reports.sort(key=lambda v: v.current_a, reverse=True)
+
+    # embedded potential + per-layer current density @ I_test; chain
+    # cells have no sheet faces in the model, so keep them out of the
+    # face computation and overlay their true 1D link density instead
+    V3 = np.full((L, ny, nx), np.nan)
+    V3[stack.masks] = Vflat[:n_grid].reshape(L, ny, nx)[stack.masks] * s
+    sheet = stack.masks if stack.chain is None \
+        else stack.masks & ~stack.chain
+    J3 = np.stack([
+        _face_current_density(
+            np.nan_to_num(V3[li]), sheet[li], sigmas[li],
+            h_m, problem.layers[li].thickness_nm * 1e-9,
+            sig2d=_sigma_2d(stack, li, sigmas[li], sigma_buildup),
+            rho=problem.rho_ohm_m)
+        for li in range(L)
+    ])
+    overlay_chain_density(stack, problem.rho_ohm_m, V3, J3)
+    return Pe, Ie, P_layers, P_vias, Parea, via_reports, V3, J3
+
+
 def _conductance_params(problem: Problem, stack: RasterStack,
                         freq_hz: float):
     """Effective (possibly AC) sheet conductances per layer, Rs ratios,
@@ -711,18 +952,9 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray,
     t0 = time.perf_counter()
     s = i_test * volts_per_amp              # unit-drive volts -> volts @ I_test
 
-    # per-edge power @ I_test; distribute in-plane power to endpoint cells
-    Pe = edges.w * ((Vflat[edges.a] - Vflat[edges.b]) * s) ** 2
-    inplane = edges.via_index < 0
-    Pflat = np.zeros(Vflat.size)
-    np.add.at(Pflat, edges.a[inplane], 0.5 * Pe[inplane])
-    np.add.at(Pflat, edges.b[inplane], 0.5 * Pe[inplane])
-    Parea = Pflat.reshape(L, ny, nx) / (h_m * h_m)
-    Parea[~stack.masks] = np.nan
-    plane = ny * nx
-    P_layers = [float(Pflat[li * plane:(li + 1) * plane].sum())
-                for li in range(L)]
-    P_vias = float(Pe[~inplane].sum())
+    Pe, Ie, P_layers, P_vias, Parea, via_reports, V3, J3 = \
+        _postprocess_fields(problem, stack, edges, Vflat, s, sigmas,
+                            sigma_buildup)
     P_total = i_test ** 2 * R
     balance = abs((sum(P_layers) + P_vias) - P_total) / max(P_total, 1e-300)
     if not np.isfinite(balance) or balance > 1e-3:
@@ -733,22 +965,6 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray,
             f"different grid size."
         )
 
-    # via reports: max segment current + total power per via
-    Ie = edges.w * (Vflat[edges.a] - Vflat[edges.b])   # amps at unit drive
-    via_reports = []
-    if problem.vias:
-        vidx = edges.via_index
-        for vi in np.unique(vidx[vidx >= 0]):
-            sel = vidx == vi
-            via = problem.vias[vi]
-            via_reports.append(ViaReport(
-                x_mm=via.x * 1e-6, y_mm=via.y * 1e-6, kind=via.kind,
-                drill_mm=via.drill_nm * 1e-6,
-                current_a=float(np.abs(Ie[sel]).max()) * s,
-                power_w=float(Pe[sel].sum()),
-            ))
-        via_reports.sort(key=lambda v: v.current_a, reverse=True)
-
     # per-injection-area currents
     part_currents1 = _part_currents(
         parts1 or [], Ie, edges, e1.ravel(), s, i_test,
@@ -756,23 +972,6 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray,
     part_currents2 = _part_currents(
         parts2 or [], Ie, edges, e2.ravel(), s, i_test,
         contact_model, int(e2.sum()))
-
-    # embedded potential + per-layer current density @ I_test; chain
-    # cells have no sheet faces in the model, so keep them out of the
-    # face computation and overlay their true 1D link density instead
-    V3 = np.full((L, ny, nx), np.nan)
-    V3[stack.masks] = Vflat.reshape(L, ny, nx)[stack.masks] * s
-    sheet = stack.masks if stack.chain is None \
-        else stack.masks & ~stack.chain
-    J3 = np.stack([
-        _face_current_density(
-            np.nan_to_num(V3[li]), sheet[li], sigmas[li],
-            h_m, problem.layers[li].thickness_nm * 1e-9,
-            sig2d=_sigma_2d(stack, li, sigmas[li], sigma_buildup),
-            rho=problem.rho_ohm_m)
-        for li in range(L)
-    ])
-    overlay_chain_density(stack, problem.rho_ohm_m, V3, J3)
     timings["postprocess_s"] = time.perf_counter() - t0
 
     return Result(
@@ -790,3 +989,551 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray,
         rs_ratios=rs_ratios,
         timings=timings,
     )
+
+
+# --- PDN mode ---------------------------------------------------------------
+#
+# N supplies + M loads instead of one driven terminal pair, solved in
+# ABSOLUTE volts (no unit-drive rescale). Each supply is a Thevenin
+# source: a virtual node held Dirichlet at v_oc, attached to its
+# contact cells through 1/(r_out * n_cells) each (sum = 1/r_out) -
+# because Dirichlet nodes are eliminated, virtual nodes never enter the
+# matrix, they only shift the diagonal/RHS of their contact cells and
+# the system stays SPD. r_out <= PDN_R_OUT_EPS degrades to a direct
+# Dirichlet contact (the exact limit). Each load draws its prescribed
+# current with uniform orthogonal injection (-I/n per cell), the same
+# semantics as the classic "uniform" contact model. Supply currents are
+# OUTCOMES (the Thevenin split); KCL makes them sum to the load draws.
+
+def _label_terminals(terminals: list) -> None:
+    """Assign S1/L1-style display tags to unlabeled terminals (in
+    definition order, per role)."""
+    ns = nl = 0
+    for t in terminals:
+        if t.role == "supply":
+            ns += 1
+            if not t.label:
+                t.label = f"S{ns}"
+        else:
+            nl += 1
+            if not t.label:
+                t.label = f"L{nl}"
+
+
+def _validate_terminals(terminals: list) -> None:
+    for t in terminals:
+        if t.role not in ("supply", "load"):
+            raise ElectrodeError(
+                f"Terminal '{t.label}': unknown role '{t.role}' "
+                f"(expected 'supply' or 'load')."
+            )
+        if t.role == "load" and t.i_draw_a < 0:
+            raise ElectrodeError(
+                f"Load '{t.label}': i_draw_a must be >= 0 "
+                f"(got {t.i_draw_a:g})."
+            )
+        if t.role == "supply" and t.r_out_ohm < 0:
+            raise ElectrodeError(
+                f"Supply '{t.label}': r_out_ohm must be >= 0 "
+                f"(got {t.r_out_ohm:g})."
+            )
+    if not any(t.role == "supply" for t in terminals):
+        raise ElectrodeError("PDN mode needs at least one supply terminal.")
+    if not any(t.role == "load" for t in terminals):
+        raise ElectrodeError("PDN mode needs at least one load terminal.")
+
+
+@dataclass
+class _Attach:
+    """How one supply is wired into the extended graph."""
+    v_oc: float
+    r_out: float
+    nodes: np.ndarray             # contact node ids (base space)
+    ideal: bool                   # r_out below eps: direct Dirichlet
+    e0: int = 0                   # its attachment edges in edges_ext
+    e1: int = 0                   # (empty slice for ideal supplies)
+
+
+def _pdn_attach(terminals: list, term_nodes: list, state: np.ndarray,
+                edges: Edges, v_nominal: float):
+    """Extend the copper graph with the PDN boundary conditions. state
+    is uint8 over the base node space (1 = copper); term_nodes carries
+    each terminal's contact node ids in that space (uniform grid: flat
+    cell ids; adaptive: leaf ids - contact cells are pinned fine there,
+    so nodes and cells are 1:1 and per-node injection equals per-cell).
+
+    A BONDED terminal shorts all its contact cells into one super-node:
+    `merge` maps every node id to its representative (identity outside
+    bonded terminals; None when no terminal is bonded). The caller
+    solves on merge-relabeled edges (member cells leave the system,
+    state 0) and afterwards scatters the potentials back with
+    Vflat = Vflat[merge], so all extraction runs on the ORIGINAL edge
+    endpoints where internal member-member edges carry exactly zero.
+
+    Returns (state_ext, dirichlet_v, inj, edges_ext, attaches, merge)
+    with virtual supply nodes appended after state.size."""
+    n_base = state.size
+    n_virt = sum(1 for t, nd in zip(terminals, term_nodes)
+                 if t.role == "supply" and len(nd)
+                 and t.r_out_ohm > config.PDN_R_OUT_EPS)
+    n_ext = n_base + n_virt
+    state_ext = np.zeros(n_ext, dtype=np.uint8)
+    state_ext[:n_base] = state
+    dirichlet_v = np.zeros(n_ext)
+    inj = np.zeros(n_ext)
+    merge = None
+
+    def bond(nodes):
+        """Short the cells to nodes[0]; members leave the system."""
+        nonlocal merge
+        if merge is None:
+            merge = np.arange(n_ext, dtype=np.int64)
+        rep = int(nodes[0])
+        merge[nodes] = rep
+        state_ext[nodes] = 0
+        return rep
+
+    aa = [edges.a]
+    bb = [edges.b]
+    ww = [edges.w]
+    vv = [edges.via_index]
+    attaches: list = []
+    nv = 0
+    e_next = len(edges.a)
+    for t, nodes in zip(terminals, term_nodes):
+        if t.role != "supply":
+            attaches.append(None)
+            if len(nodes):
+                if t.bonded:
+                    rep = bond(nodes)
+                    state_ext[rep] = 1
+                    inj[rep] -= t.i_draw_a      # whole draw at the lug
+                else:
+                    inj[nodes] -= t.i_draw_a / len(nodes)
+            continue
+        v = v_nominal if t.v_oc is None else t.v_oc
+        at = _Attach(v_oc=v, r_out=t.r_out_ohm, nodes=nodes,
+                     ideal=t.r_out_ohm <= config.PDN_R_OUT_EPS)
+        attaches.append(at)
+        if len(nodes) == 0:
+            continue                        # restricted away: reports 0 A
+        if t.bonded:
+            rep = bond(nodes)
+            if at.ideal:
+                state_ext[rep] = 2
+                dirichlet_v[rep] = v
+            else:
+                state_ext[rep] = 1
+                vid = n_base + nv
+                nv += 1
+                state_ext[vid] = 2
+                dirichlet_v[vid] = v
+                # the whole r_out in series with the equipotential lug
+                aa.append(np.array([vid], dtype=np.int64))
+                bb.append(np.array([rep], dtype=np.int64))
+                ww.append(np.array([1.0 / t.r_out_ohm]))
+                vv.append(np.array([PDN_EDGE], dtype=np.int32))
+                at.e0, at.e1 = e_next, e_next + 1
+                e_next += 1
+        elif at.ideal:
+            state_ext[nodes] = 2
+            dirichlet_v[nodes] = v
+        else:
+            vid = n_base + nv
+            nv += 1
+            state_ext[vid] = 2
+            dirichlet_v[vid] = v
+            k = len(nodes)
+            # oriented virtual -> cell so Ie = w * (v_oc - V_cell) is
+            # the delivered current, positive out of the supply
+            aa.append(np.full(k, vid, dtype=np.int64))
+            bb.append(nodes.astype(np.int64))
+            ww.append(np.full(k, 1.0 / (t.r_out_ohm * k)))
+            vv.append(np.full(k, PDN_EDGE, dtype=np.int32))
+            at.e0, at.e1 = e_next, e_next + k
+            e_next += k
+    edges_ext = Edges(a=np.concatenate(aa), b=np.concatenate(bb),
+                      w=np.concatenate(ww), via_index=np.concatenate(vv),
+                      dead_barrels=edges.dead_barrels)
+    return state_ext, dirichlet_v, inj, edges_ext, attaches, merge
+
+
+def _pdn_solve_edges(edges_ext: Edges, merge) -> Edges:
+    """The edge set the linear system is assembled from: bonded
+    terminals' member cells relabeled to their representative (edge
+    ORDER and COUNT are preserved - only endpoints move; internal
+    edges become self-loops, which cancel exactly in the COO
+    assembly). Identity when nothing is bonded."""
+    if merge is None:
+        return edges_ext
+    return Edges(a=merge[edges_ext.a], b=merge[edges_ext.b],
+                 w=edges_ext.w, via_index=edges_ext.via_index,
+                 dead_barrels=edges_ext.dead_barrels)
+
+
+def _pdn_extract(terminals: list, term_nodes: list, attaches: list,
+                 Vflat: np.ndarray, Ie: np.ndarray, edges_ext: Edges,
+                 term_part_nodes: list):
+    """Per-supply delivered currents and per-load voltages/powers from
+    the solved potentials, shared by the uniform-grid and adaptive PDN
+    paths (Ie must satisfy KCL - the corrected currents on the adaptive
+    path; bonded terminals' Vflat already scattered back, so their
+    internal edges carry exactly zero). Returns (supplies, loads)."""
+    n_ext = Vflat.size
+    copper = edges_ext.via_index != PDN_EDGE
+
+    def part_flux_out(pn):
+        """Oriented copper-edge flux out of a part's cells: the part's
+        boundary current (same-terminal internal edges carry zero for
+        Dirichlet and bonded contacts; attachment edges excluded so a
+        bonded supply's lug edge is not double-counted)."""
+        pm = np.zeros(n_ext, dtype=bool)
+        pm[pn] = True
+        return float(Ie[pm[edges_ext.a] & copper].sum()
+                     - Ie[pm[edges_ext.b] & copper].sum())
+
+    supplies: list = []
+    loads: list = []
+    for t, nodes, at, pnodes in zip(terminals, term_nodes, attaches,
+                                    term_part_nodes):
+        if t.role == "supply":
+            if len(at.nodes) == 0:
+                supplies.append(SupplyReport(
+                    label=t.label, v_oc=at.v_oc, r_out_ohm=at.r_out,
+                    i_a=0.0, v_contact=at.v_oc, p_internal_w=0.0,
+                    part_currents=[(pl, 0.0) for pl, _ in pnodes],
+                    v_eff=at.v_oc, component=t.component,
+                    comment=t.comment))
+                continue
+            if at.ideal:
+                # exact discrete flux out of the Dirichlet contact
+                # (edges inside the contact cancel a-side vs b-side)
+                member = np.zeros(n_ext, dtype=bool)
+                member[at.nodes] = True
+                i_a = float(Ie[member[edges_ext.a] & copper].sum()
+                            - Ie[member[edges_ext.b] & copper].sum())
+                v_contact = at.v_oc
+                v_eff = at.v_oc
+                p_int = 0.0
+                pcs = [(pl, part_flux_out(pn)) for pl, pn in pnodes]
+            else:
+                sl = slice(at.e0, at.e1)
+                i_a = float(Ie[sl].sum())
+                v_contact = float(Vflat[at.nodes].mean())
+                dv = Vflat[edges_ext.a[sl]] - Vflat[edges_ext.b[sl]]
+                p_int = float((edges_ext.w[sl] * dv * dv).sum())
+                if t.bonded:
+                    # one lug edge carries the total; the per-part
+                    # split is the copper boundary flux (an outcome)
+                    v_eff = v_contact       # equipotential lug
+                    pcs = [(pl, part_flux_out(pn)) for pl, pn in pnodes]
+                else:
+                    cells = edges_ext.b[sl]
+                    # the potential the delivered power actually sees:
+                    # per-cell currents weight their cell potentials
+                    v_eff = (float((Ie[sl] * Vflat[cells]).sum() / i_a)
+                             if abs(i_a) > 1e-300 else v_contact)
+                    pcs = []
+                    for pl, pn in pnodes:
+                        pm = np.zeros(n_ext, dtype=bool)
+                        pm[pn] = True
+                        pcs.append((pl, float(Ie[sl][pm[cells]].sum())))
+            supplies.append(SupplyReport(
+                label=t.label, v_oc=at.v_oc, r_out_ohm=at.r_out,
+                i_a=i_a, v_contact=v_contact, p_internal_w=p_int,
+                part_currents=pcs, v_eff=v_eff,
+                component=t.component, comment=t.comment))
+        else:
+            v = Vflat[nodes]
+            n_k = len(nodes)
+            if t.bonded:
+                # split by the network through the external bond
+                pcs = [(pl, -part_flux_out(pn)) for pl, pn in pnodes]
+            else:
+                pcs = [(pl, t.i_draw_a * len(pn) / max(n_k, 1))
+                       for pl, pn in pnodes]
+            loads.append(LoadReport(
+                label=t.label, i_a=t.i_draw_a,
+                v_mean=float(v.mean()), v_min=float(v.min()),
+                p_w=t.i_draw_a * float(v.mean()), part_currents=pcs,
+                component=t.component, comment=t.comment))
+    return supplies, loads
+
+
+def _pdn_balance(supplies: list, loads: list, P_layers: list,
+                 P_vias: float) -> tuple[float, float, float, float, float]:
+    """Generalized power balance (Tellegen): source power = copper
+    dissipation + R_out dissipation + load power, exactly, independent
+    of the voltage reference (supply and load currents sum equal).
+    Returns (balance_rel, mismatch_rel, i_sup, i_loads, p_loads_rout)
+    or raises SolverError; a zero-power probe run (no draws, equal
+    v_oc) has nothing to balance and reports 0."""
+    i_loads = sum(l.i_a for l in loads)
+    i_sup = sum(s_.i_a for s_ in supplies)
+    p_rout = sum(s_.p_internal_w for s_ in supplies)
+    p_loads = sum(l.p_w for l in loads)
+    p_src = sum(s_.v_oc * s_.i_a for s_ in supplies)
+    # gross scale, not |net|: with circulating currents between supplies
+    # the net source power nearly cancels while watts really flow - the
+    # residual must be judged against what flows, or the check trips on
+    # pure floating-point cancellation
+    p_gross = sum(abs(s_.v_oc * s_.i_a) for s_ in supplies)
+    p_sink = sum(P_layers) + P_vias + p_rout + p_loads
+    vocs = [s_.v_oc for s_ in supplies]
+    if i_loads == 0.0 and max(vocs) == min(vocs):
+        balance = 0.0
+    else:
+        balance = abs(p_src - p_sink) / max(p_gross, 1e-300)
+        if not np.isfinite(balance) or balance > 1e-3:
+            raise SolverError(
+                f"Inconsistent PDN solve: power-balance error "
+                f"{balance:.2e} (sources {p_src:.6g} W vs copper + "
+                f"R_out + loads {p_sink:.6g} W). The result is not "
+                f"trustworthy - try a different grid size."
+            )
+    i_scale = max(abs(i_loads),
+                  max((abs(s_.i_a) for s_ in supplies), default=0.0),
+                  1e-300)
+    mismatch = abs(i_sup - i_loads) / i_scale
+    return balance, mismatch, i_sup, i_loads, p_loads
+
+
+def _pdn_pairs(terminals: list, term_nodes: list, attaches: list,
+               merge, state_base: np.ndarray, edges: Edges,
+               supplies: list, loads: list, make_solver) -> list:
+    """Effective copper resistance between every supply and every load,
+    plus the proportional-sharing loss allocation (see PairReport).
+
+    The pair network is the COPPER alone - attachment resistances and
+    Thevenin sources stripped. Contact patterns mirror the solve's
+    models: loads and resistive supplies inject uniformly per cell,
+    bonded terminals and ideal supplies are equipotential super-nodes.
+    R_ij = (p_i - p_j)^T L_g^{-1} (p_i - p_j) costs ONE extra solve per
+    terminal (same factorization; one node per connected component is
+    grounded - the balanced quadratic form is ground-independent).
+    Pairs without a common copper component get r_ohm None and no
+    allocation; the allocation splits each component's copper loss by
+    f_ij = I_i * I_j / I_loads_of_that_component, which sums exactly
+    to the total copper dissipation.
+
+    make_solver(state_g, dirichlet_v, edges_pm, pmerge) -> solve(inj)
+    -> V lets the adaptive path plug in its deferred-correction loop
+    (pmerge is the pair system's own node-merge map, or None)."""
+    n_base = state_base.size
+    idx0 = np.arange(n_base, dtype=np.int64)
+    # pair-system merge: the solve's bonded lugs plus every ideal
+    # supply contact (a Dirichlet region is equipotential too)
+    pmerge = idx0.copy() if merge is None else merge[:n_base].copy()
+    for t, at, nodes in zip(terminals, attaches, term_nodes):
+        if (t.role == "supply" and at is not None and at.ideal
+                and not t.bonded and len(nodes)):
+            pmerge[nodes] = pmerge[int(nodes[0])]
+    if bool((pmerge == idx0).all()):
+        pmerge = None
+
+    state_g = state_base.copy()
+    if pmerge is not None:
+        state_g[pmerge != idx0] = 0
+        edges_pm = Edges(a=pmerge[edges.a], b=pmerge[edges.b],
+                         w=edges.w, via_index=edges.via_index,
+                         dead_barrels=edges.dead_barrels)
+    else:
+        edges_pm = edges
+
+    pats = []
+    for t, at, nodes in zip(terminals, attaches, term_nodes):
+        if len(nodes) == 0:
+            pats.append(None)
+            continue
+        p = np.zeros(n_base)
+        if (t.bonded or (t.role == "supply" and at is not None
+                         and at.ideal)):
+            rep = (int(pmerge[nodes[0]]) if pmerge is not None
+                   else int(nodes[0]))
+            p[rep] = 1.0
+        else:
+            p[nodes] = 1.0 / len(nodes)
+        pats.append(p)
+
+    graph = sparse.coo_matrix(
+        (np.ones(len(edges_pm.a)), (edges_pm.a, edges_pm.b)),
+        shape=(n_base, n_base))
+    _, labels = csgraph.connected_components(graph, directed=False)
+    tcomp = []
+    for p in pats:
+        if p is None:
+            tcomp.append(None)
+            continue
+        cs = set(labels[np.flatnonzero(p)].tolist())
+        # a terminal spread over several components has no single
+        # pair resistance - its pairs report "no path"
+        tcomp.append(cs.pop() if len(cs) == 1 else None)
+
+    dv = np.zeros(n_base)
+    grounds: set = set()
+    for c, p in zip(tcomp, pats):
+        if c is not None and c not in grounds:
+            state_g[int(np.flatnonzero(p)[0])] = 2
+            grounds.add(c)
+    if not grounds:
+        return []
+
+    progress.stage("source-sink pair resistances ...")
+    solve = make_solver(state_g, dv, edges_pm, pmerge)
+    us = [solve(p) if p is not None and c is not None else None
+          for p, c in zip(pats, tcomp)]
+
+    # per-component load current: proportional sharing never crosses a
+    # copper gap (nothing flows between disconnected sheets)
+    comp_load_i: dict = {}
+    lidx = -1
+    for j, tj in enumerate(terminals):
+        if tj.role != "load":
+            continue
+        lidx += 1
+        if tcomp[j] is not None:
+            comp_load_i[tcomp[j]] = (comp_load_i.get(tcomp[j], 0.0)
+                                     + loads[lidx].i_a)
+
+    rows: list = []
+    sidx = -1
+    for i, ti in enumerate(terminals):
+        if ti.role != "supply":
+            continue
+        sidx += 1
+        s = supplies[sidx]
+        lidx = -1
+        for j, tj in enumerate(terminals):
+            if tj.role != "load":
+                continue
+            lidx += 1
+            ld = loads[lidx]
+            same = tcomp[i] is not None and tcomp[i] == tcomp[j]
+            r = None
+            f = 0.0
+            if same:
+                r = float(pats[i] @ us[i] + pats[j] @ us[j]
+                          - 2.0 * (pats[i] @ us[j]))
+                i_tot = comp_load_i.get(tcomp[i], 0.0)
+                if i_tot > 0.0:
+                    f = s.i_a * ld.i_a / i_tot
+            rows.append(PairReport(
+                supply=s.label, load=ld.label, r_ohm=r,
+                i_share_a=f, p_w=f * (s.v_eff - ld.v_mean)))
+    return rows
+
+
+def run_solve_pdn(problem: Problem, stack: RasterStack, term_masks: list,
+                  term_parts: list, freq_hz: float = 0.0,
+                  v_nominal: float | None = None) -> Result:
+    """PDN solve on the uniform grid (adaptive dispatch on top, like
+    run_solve). term_masks/term_parts come from raster.terminal_masks /
+    terminal_partition, aligned with problem.terminals. Contact models
+    are fixed: Thevenin supplies, uniform-injection loads."""
+    if v_nominal is None:
+        v_nominal = config.PDN_V_NOMINAL
+    terminals = problem.terminals
+    _label_terminals(terminals)
+    _validate_terminals(terminals)
+    if config.ADAPTIVE_CELLS:
+        from . import adaptive
+        return adaptive.run_solve_adaptive_pdn(problem, stack, term_masks,
+                                               term_parts, freq_hz,
+                                               v_nominal)
+
+    timings = {}
+    sigmas, rs_ratios, via_factor, sigma_buildup = \
+        _conductance_params(problem, stack, freq_hz)
+
+    t0 = time.perf_counter()
+    edges = build_edges(stack, problem, sigmas, via_factor, sigma_buildup)
+    changed, _n_kept = connected_restrict_multi(stack, term_masks,
+                                                terminals, edges)
+    if changed:
+        edges = build_edges(stack, problem, sigmas, via_factor,
+                            sigma_buildup)
+    if edges.dead_barrels:
+        print(f"warning: {edges.dead_barrels} via/pad barrel(s) found fill "
+              f"copper on fewer than 2 layers and carry no current (pad "
+              f"copper is not modeled; a finer grid may pick up thermal "
+              f"spokes)")
+    if stack.buildup is not None:
+        stack.buildup &= stack.masks
+    if stack.chain is not None:
+        stack.chain &= stack.masks
+    for t, parts in zip(terminals, term_parts):
+        for label, m in parts:
+            had = bool(m.any())
+            m &= stack.masks
+            if had and not m.any():
+                print(f"warning: contact part '{label}' of {t.role} "
+                      f"'{t.label}' only touches disconnected copper - "
+                      f"it carries no current")
+    timings["edges_s"] = time.perf_counter() - t0
+
+    t0 = time.perf_counter()
+    state = np.zeros(stack.masks.size, dtype=np.uint8)
+    state[stack.masks.ravel()] = 1
+    term_nodes = [np.flatnonzero(m.ravel()) for m in term_masks]
+    state_base = state.copy()       # copper-only state for _pdn_pairs
+    state, dirichlet_v, inj, edges_ext, attaches, merge = _pdn_attach(
+        terminals, term_nodes, state, edges, v_nominal)
+    A, rhs, _ = _assemble(state, _pdn_solve_edges(edges_ext, merge),
+                          inj, dirichlet_v)
+    x, info = solve_system(A, rhs)
+    Vflat = np.where(state >= 2, dirichlet_v, 0.0)
+    Vflat[state == 1] = x
+    if merge is not None:
+        Vflat = Vflat[merge]        # bonded members read their lug
+    timings["solve_s"] = time.perf_counter() - t0
+
+    t0 = time.perf_counter()
+    Pe, Ie, P_layers, P_vias, Parea, via_reports, V3, J3 = \
+        _postprocess_fields(problem, stack, edges_ext, Vflat, 1.0, sigmas,
+                            sigma_buildup)
+    term_part_nodes = [
+        [(pl, np.flatnonzero(m.ravel())) for pl, m in parts]
+        for parts in term_parts]
+    supplies, loads = _pdn_extract(terminals, term_nodes, attaches, Vflat,
+                                   Ie, edges_ext, term_part_nodes)
+    balance, mismatch, i_sup, i_loads, p_loads = _pdn_balance(
+        supplies, loads, P_layers, P_vias)
+    timings["postprocess_s"] = time.perf_counter() - t0
+
+    t0 = time.perf_counter()
+
+    def _pair_solver(state_g, dv, edges_pm, pmerge):
+        A2, rhs0p, _ = _assemble(state_g, edges_pm, None, dv)
+        ps2 = PreparedSolver(A2)
+        freeg = state_g == 1
+
+        def slv(inj_p):
+            x2, _ = ps2.solve(rhs0p + inj_p[freeg])
+            V = np.where(state_g >= 2, dv, 0.0)
+            V[freeg] = x2
+            return V
+        return slv
+
+    pairs = _pdn_pairs(terminals, term_nodes, attaches, merge,
+                       state_base, edges, supplies, loads, _pair_solver)
+    timings["pairs_s"] = time.perf_counter() - t0
+
+    return Result(
+        R_ohm=float("nan"), i_test=i_loads, V=V3, Jmag=J3, Parea=Parea,
+        layer_names=list(stack.layer_names),
+        P_total=float(sum(P_layers) + P_vias),
+        P_layers=P_layers, P_vias=P_vias,
+        power_balance_rel=balance, via_reports=via_reports,
+        I1_a=i_sup, I2_a=i_loads, mismatch_rel=mismatch,
+        n_free=info.n_unknowns, solve_info=info,
+        contact_model="pdn",
+        freq_hz=freq_hz,
+        skin_depth_um=(skin.skin_depth_m(freq_hz, problem.rho_ohm_m) * 1e6
+                       if freq_hz > 0 else None),
+        rs_ratios=rs_ratios,
+        timings=timings,
+        mode="pdn", supplies=supplies, loads=loads,
+        P_loads=p_loads,
+        P_supply_internal=sum(s_.p_internal_w for s_ in supplies),
+        v_nominal=v_nominal, pairs=pairs,
+    )