diff --git a/README.md b/README.md index af5b5eb..a8d43cb 100644 --- a/README.md +++ b/README.md @@ -137,6 +137,17 @@ SWIG API. Requires KiCad **10.0.1+**. AMG-preconditioned CG (pyamg) above — Jacobi-CG if pyamg is missing. Discretization error typically ≲ 2 % at defaults — halve the cell size and compare to judge convergence. +- **Adaptive cells** (dialog checkbox, off by default; `ADAPTIVE_CELLS`): + solves on a 2:1-balanced quadtree — fine cells at copper boundaries, + electrodes, traces, via mouths and buildup, blocks up to + `ADAPTIVE_MAX_CELL_UM` (2 mm) in plane interiors (`ADAPTIVE_GUARD` + sets the clearance a block needs to grow). The **minimum element size + is the grid cell size itself** (auto / dialog / `CELL_UM_OVERRIDE`); + the uniform limit reproduces the normal grid exactly. Large + speed/memory wins on big pours; the coarse–fine interfaces carry a + first-order flux error that biases R **low by ~0.5–2 %** depending on + geometry (worst on narrow strips, mild on large planes). All fields + are expanded back to the fine grid for the maps and reports. ## Offline / development diff --git a/fill_resistance/adaptive.py b/fill_resistance/adaptive.py new file mode 100644 index 0000000..fc116b6 --- /dev/null +++ b/fill_resistance/adaptive.py @@ -0,0 +1,348 @@ +"""Adaptive-grid solve path (phase 2): maps the fully rasterized problem +onto per-layer balanced quadtree leaf graphs (quadtree.py), solves with +the production assembly/AMG, and expands every field back to the fine +grid, so plots, reports and dumps are unchanged. + +Enabled via config.ADAPTIVE_CELLS (dialog checkbox "adaptive cells"). +Every fine cell that carries anything non-uniform - electrodes, 1D +trace-chain cells, solder buildup, via-mouth thickness scaling - is +pinned at the fine size (keep_fine), so all coarser leaves have the +plain layer conductance and the leaf system reduces EXACTLY to the +production system wherever the grid is fine. The minimum element size +is therefore the grid cell size itself; ADAPTIVE_MAX_CELL_UM caps the +coarsest leaf. + +ACCURACY: coarse-fine interfaces carry a first-order two-point-flux +error (centers of different-size neighbors are laterally offset), which +biases R LOW by ~0.5-2% depending on geometry - worst where transition +rings span much of the current path (narrow strips), mild on large +pours. Symmetric fine pairs under a coarse face cancel pairwise; the +residue comes from unpaired larger-neighbor faces. Gradient-corrected +interface fluxes (phase 4) are the known cure if tighter accuracy per +leaf is ever needed. +""" +from __future__ import annotations + +import time + +import numpy as np +from scipy import sparse +from scipy.sparse import csgraph + +from . import config, quadtree, skin +from . import solver as sv +from .errors import ConnectivityError +from .geometry import Problem +from .raster import RasterStack + + +def _max_block(h_nm: float) -> int: + mb = 1 + while mb * 2 * h_nm <= config.ADAPTIVE_MAX_CELL_UM * 1000.0: + mb *= 2 + return mb + + +def _nodes_of_cells(grids, offs, li: int, cells2d: np.ndarray) -> np.ndarray: + """Node ids of the (copper) fine cells selected by a 2D bool mask.""" + ids = grids[li].id_grid[cells2d] + ids = ids[ids >= 0].astype(np.int64) + return offs[li] + np.unique(ids) + + +def run_solve_adaptive(problem: Problem, stack: RasterStack, + e1: np.ndarray, e2: np.ndarray, i_test: float, + freq_hz: float, contact_model: str, + parts1: list | None, + parts2: list | None) -> sv.Result: + timings = {} + L, ny, nx = stack.masks.shape + h_m = stack.h_nm * 1e-9 + plane = ny * nx + + sigmas, rs_ratios, via_factor, sigma_buildup = \ + sv._conductance_params(problem, stack, freq_hz) + + # --- leaves per layer ------------------------------------------------- + t0 = time.perf_counter() + keep = e1 | e2 + if stack.chain is not None: + keep |= stack.chain + if stack.buildup is not None: + keep |= stack.buildup + if stack.thick_scale is not None: + keep |= stack.thick_scale != 1.0 + mb = _max_block(stack.h_nm) + grids = [quadtree.build_leaves(stack.masks[li], keep_fine=keep[li], + max_block=mb, + guard=config.ADAPTIVE_GUARD) + for li in range(L)] + offs = np.zeros(L + 1, dtype=np.int64) + for li in range(L): + offs[li + 1] = offs[li] + grids[li].n + N = int(offs[-1]) + n_cells = int(stack.masks.sum()) + print(f"adaptive grid: {N} leaves for {n_cells} copper cells " + f"({n_cells / max(N, 1):.1f}x, max leaf " + f"{max(int(g.size.max()) if g.n else 1 for g in grids)} cells)") + + # --- edges: in-plane faces, 1D chain links, barrels ------------------- + aa, bb, ww, vv = [], [], [], [] + sig_leaves, teq_leaves = [], [] + for li in range(L): + g_ = grids[li] + sig_leaf = np.full(g_.n, sigmas[li]) + t_m = problem.layers[li].thickness_nm * 1e-9 + s2d = sv._sigma_2d(stack, li, sigmas[li], sigma_buildup) + fine = g_.size == 1 + if s2d is not None and fine.any(): + sig_leaf[fine] = s2d[g_.y0[fine], g_.x0[fine]] + # J reference thickness: same convention as the uniform grid + teq_leaves.append(sig_leaf * problem.rho_ohm_m if s2d is not None + else np.full(g_.n, t_m)) + sig_leaves.append(sig_leaf) + chainleaf = np.zeros(g_.n, dtype=bool) + if stack.chain is not None and fine.any(): + chainleaf[fine] = stack.chain[li][g_.y0[fine], g_.x0[fine]] + ia, ib, wl, ax = quadtree.leaf_faces(g_) + ok = ~(chainleaf[ia] | chainleaf[ib]) + ia, ib, wl, ax = ia[ok], ib[ok], wl[ok], ax[ok] + gcond = wl / (g_.size[ia] / (2.0 * sig_leaf[ia]) + + g_.size[ib] / (2.0 * sig_leaf[ib])) + aa.append(offs[li] + ia) + bb.append(offs[li] + ib) + ww.append(gcond) + vv.append(np.full(len(ia), -1, dtype=np.int32)) + + if stack.chain_edges is not None and len(stack.chain_edges[0]): + ca, cb, cg, cl = stack.chain_edges + alive = stack.masks.ravel()[ca] & stack.masks.ravel()[cb] + if alive.any(): + na = np.empty(len(ca), dtype=np.int64) + nb = np.empty(len(ca), dtype=np.int64) + for flat, out in ((ca, na), (cb, nb)): + li_ = flat // plane + rem = flat - li_ * plane + for l in range(L): + m = li_ == l + if m.any(): + ids = grids[l].id_grid[rem[m] // nx, rem[m] % nx] + out[m] = np.where(ids >= 0, offs[l] + ids, -1) + alive &= (na >= 0) & (nb >= 0) + fac = np.array([sigmas[l] * problem.rho_ohm_m + / (problem.layers[l].thickness_nm * 1e-9) + for l in range(L)]) + aa.append(na[alive]) + bb.append(nb[alive]) + ww.append((cg * fac[cl])[alive]) + vv.append(np.full(int(alive.sum()), -1, dtype=np.int32)) + + links, dead_barrels = sv._barrel_links(stack, problem) + for vi, la, ia_, ja_, lb, ib_, jb_, r_dc in links: + na = offs[la] + grids[la].id_grid[ia_, ja_] + nb = offs[lb] + grids[lb].id_grid[ib_, jb_] + aa.append(np.array([na], dtype=np.int64)) + bb.append(np.array([nb], dtype=np.int64)) + ww.append(np.array([1.0 / (r_dc * via_factor)])) + vv.append(np.array([vi], dtype=np.int32)) + if dead_barrels: + print(f"warning: {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 not aa: + raise ConnectivityError("No copper found on the selected layers.") + edges = sv.Edges(a=np.concatenate(aa), b=np.concatenate(bb), + w=np.concatenate(ww), via_index=np.concatenate(vv), + dead_barrels=dead_barrels) + + # --- connectivity restriction on the leaf graph ----------------------- + graph = sparse.coo_matrix( + (np.ones(len(edges.a)), (edges.a, edges.b)), shape=(N, N)) + _, labels = csgraph.connected_components(graph, directed=False) + e1n = np.zeros(N, dtype=bool) + e2n = np.zeros(N, dtype=bool) + for li in range(L): + e1n[_nodes_of_cells(grids, offs, li, e1[li])] = True + e2n[_nodes_of_cells(grids, offs, li, e2[li])] = True + common = np.intersect1d(np.unique(labels[e1n]), np.unique(labels[e2n])) + if len(common) == 0: + raise ConnectivityError( + "The two terminals are not connected by the selected fill " + "layers (not even through vias). Check the layer selection and " + "that the fills are up to date." + ) + if len(common) > 1 and contact_model != "equipotential": + raise ConnectivityError( + f"The selected fills form {len(common)} disconnected copper " + f"groups that each touch both terminals. The uniform-injection " + f"contact model cannot determine the current split between " + f"disconnected sheets - switch to the equipotential contact " + f"model (bonded lug), or include the layers/vias that join " + f"them." + ) + keepn = np.isin(labels, common) + if not keepn.all(): + sel = keepn[edges.a] & keepn[edges.b] + edges = sv.Edges(a=edges.a[sel], b=edges.b[sel], w=edges.w[sel], + via_index=edges.via_index[sel], + dead_barrels=dead_barrels) + for li in range(L): + ids = grids[li].id_grid + kept_cells = (ids >= 0) & keepn[offs[li] + np.maximum(ids, 0)] + stack.masks[li] &= kept_cells + e1[li] &= kept_cells + e2[li] &= kept_cells + if stack.buildup is not None: + stack.buildup &= stack.masks + if stack.chain is not None: + stack.chain &= stack.masks + e1n &= keepn + e2n &= keepn + for label, m in (parts1 or []) + (parts2 or []): + had = bool(m.any()) + m &= stack.masks + if had and not m.any(): + print(f"warning: contact part '{label}' only touches copper " + f"that is not connected to both terminals - it carries " + f"no current") + timings["edges_s"] = time.perf_counter() - t0 + + # --- solve ------------------------------------------------------------- + t0 = time.perf_counter() + state = np.zeros(N, dtype=np.uint8) + state[keepn] = 1 + if contact_model == "equipotential": + state[e1n] = 2 + state[e2n] = 3 + Vflat, R, I1, I2, mismatch, volts_per_amp, info = \ + sv._equipotential_core(state, edges) + else: + n1, n2 = int(e1n.sum()), int(e2n.sum()) + inj = np.zeros(N) + inj[e1n] = 1.0 / n1 + inj[e2n] = -1.0 / n2 + ground = int(np.flatnonzero(e2n)[0]) + state[ground] = 3 + Vflat, R, I1, I2, mismatch, volts_per_amp, info = \ + sv._uniform_core(state, inj, e1n, e2n, edges) + timings["solve_s"] = time.perf_counter() - t0 + + # --- fields on leaves, expanded to the fine grid ------------------------ + t0 = time.perf_counter() + s = i_test * volts_per_amp + + Pe = edges.w * ((Vflat[edges.a] - Vflat[edges.b]) * s) ** 2 + inplane = edges.via_index < 0 + Pnode = np.zeros(N) + np.add.at(Pnode, edges.a[inplane], 0.5 * Pe[inplane]) + np.add.at(Pnode, edges.b[inplane], 0.5 * Pe[inplane]) + P_layers = [float(Pnode[offs[li]:offs[li + 1]].sum()) for li in range(L)] + P_vias = float(Pe[~inplane].sum()) + 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: + raise sv.SolverError( + f"Inconsistent solve: R = {R:.6g} ohm with power-balance error " + f"{balance:.2e} (sum of edge powers vs I^2*R). The result is " + f"not trustworthy - try the equipotential contact model or a " + f"different grid size." + ) + + 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(sv.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) + + def part_currents(parts, e_nodes, n_total_cells): + out = [] + for label, mask3 in (parts or []): + n_part = int(mask3.sum()) + if contact_model == "uniform": + amps = i_test * n_part / max(n_total_cells, 1) + else: + pf = np.zeros(N, dtype=bool) + for li in range(L): + pf[_nodes_of_cells(grids, offs, li, mask3[li])] = True + pf &= e_nodes + amps = abs(float(Ie[pf[edges.a]].sum() + - Ie[pf[edges.b]].sum())) * s + out.append((label, amps)) + return out + + part_currents1 = part_currents(parts1, e1n, int(e1.sum())) + part_currents2 = part_currents(parts2, e2n, int(e2.sum())) + + V3 = np.full((L, ny, nx), np.nan) + J3 = np.full((L, ny, nx), np.nan) + Parea = np.full((L, ny, nx), np.nan) + for li in range(L): + g_ = grids[li] + ids = g_.id_grid + m = stack.masks[li] + Vl = Vflat[offs[li]:offs[li + 1]] + V3[li][m] = Vl[ids[m]] * s + + # per-leaf |J| from face currents at unit drive, reconstructed + # with the same series-half-cell rule (edges were filtered by + # the restriction, so recompute locally) + gio = offs[li] + Ixn = np.zeros(g_.n) + Iyn = np.zeros(g_.n) + sig_leaf = sig_leaves[li] + ia2, ib2, wl2, ax2 = quadtree.leaf_faces(g_) + chain_ok = np.ones(len(ia2), dtype=bool) + if stack.chain is not None: + fine = g_.size == 1 + cl = np.zeros(g_.n, dtype=bool) + if fine.any(): + cl[fine] = stack.chain[li][g_.y0[fine], g_.x0[fine]] + chain_ok = ~(cl[ia2] | cl[ib2]) + ia2, ib2, wl2, ax2 = (ia2[chain_ok], ib2[chain_ok], wl2[chain_ok], + ax2[chain_ok]) + keep_f = keepn[gio + ia2] & keepn[gio + ib2] + ia2, ib2, wl2, ax2 = ia2[keep_f], ib2[keep_f], wl2[keep_f], \ + ax2[keep_f] + g2 = wl2 / (g_.size[ia2] / (2.0 * sig_leaf[ia2]) + + g_.size[ib2] / (2.0 * sig_leaf[ib2])) + If = g2 * (Vflat[gio + ia2] - Vflat[gio + ib2]) + for axis, acc in ((0, Ixn), (1, Iyn)): + selx = ax2 == axis + np.add.at(acc, ia2[selx], If[selx]) + np.add.at(acc, ib2[selx], If[selx]) + span_m = g_.size.astype(float) * h_m + with np.errstate(invalid="ignore", divide="ignore"): + Jl = np.hypot(0.5 * Ixn, 0.5 * Iyn) / (span_m * teq_leaves[li]) + J3[li][m] = Jl[ids[m]] * s + + cellP = Pnode[offs[li]:offs[li + 1]] \ + / (g_.size.astype(float) ** 2 * h_m * h_m) + Parea[li][m] = cellP[ids[m]] + timings["postprocess_s"] = time.perf_counter() - t0 + + return sv.Result( + R_ohm=R, i_test=i_test, V=V3, Jmag=J3, Parea=Parea, + layer_names=list(stack.layer_names), + P_total=P_total, P_layers=P_layers, P_vias=P_vias, + power_balance_rel=balance, via_reports=via_reports, + I1_a=I1, I2_a=I2, mismatch_rel=mismatch, + n_free=info.n_unknowns, solve_info=info, + part_currents1=part_currents1, part_currents2=part_currents2, + contact_model=contact_model, + 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, + ) diff --git a/fill_resistance/config.py b/fill_resistance/config.py index 0aed78a..c47d816 100644 --- a/fill_resistance/config.py +++ b/fill_resistance/config.py @@ -55,6 +55,20 @@ ELECTRODE_POS_LAYER = "User.1" # rectangles on this layer mark V+ contact parts ELECTRODE_NEG_LAYER = "User.2" # rectangles on this layer mark V- contact parts ALWAYS_REFILL = False # refill zones even if KiCad says they are filled +# --- Adaptive grid --- +ADAPTIVE_CELLS = False # solve on a 2:1-balanced quadtree: fine at + # copper boundaries/electrodes/features, + # coarse plane interiors (dialog-toggleable). + # Slight low bias (~0.5-1% on feature-dense + # boards); large speed/memory wins on pours +ADAPTIVE_MAX_CELL_UM = 2000.0 # coarsest leaf edge length. The MINIMUM + # element size is the grid cell size itself + # (auto / dialog / CELL_UM_OVERRIDE). Rarely + # needs tuning: the guard distance already + # caps leaf growth near features +ADAPTIVE_GUARD = 4 # a leaf of size s needs >= GUARD*s cells of + # clearance to the nearest feature + # --- Solver --- CONTACT_MODEL = "uniform" # "uniform": conductor pressed on top injects # orthogonally with uniform surface density diff --git a/fill_resistance/dialog.py b/fill_resistance/dialog.py index d8652b7..e8c08c4 100644 --- a/fill_resistance/dialog.py +++ b/fill_resistance/dialog.py @@ -37,6 +37,7 @@ class Selection: extra_cu_um: float = 0.0 include_tracks: bool = True vias_capped: bool = True + adaptive: bool = False class _Dialog(QDialog): @@ -72,6 +73,12 @@ class _Dialog(QDialog): self.capped_check.setChecked(config.VIAS_CAPPED) form.addRow("Vias:", self.capped_check) + self.adaptive_check = QCheckBox( + "adaptive cells (coarsen plane interiors; faster on large " + "boards, ~0.5–2 % low bias)") + self.adaptive_check.setChecked(config.ADAPTIVE_CELLS) + form.addRow("Grid:", self.adaptive_check) + self.contact1_box = QComboBox() self.contact2_box = QComboBox() form.addRow(f"V+ ({e1_label}):", self.contact1_box) @@ -214,7 +221,8 @@ class _Dialog(QDialog): include_buildup=self.buildup_check.isChecked(), extra_cu_um=extra_cu, include_tracks=self.tracks_check.isChecked(), - vias_capped=self.capped_check.isChecked()) + vias_capped=self.capped_check.isChecked(), + adaptive=self.adaptive_check.isChecked()) def _try_accept(self) -> None: try: diff --git a/fill_resistance/main.py b/fill_resistance/main.py index ad395e0..3cd63ac 100644 --- a/fill_resistance/main.py +++ b/fill_resistance/main.py @@ -86,6 +86,7 @@ def main() -> None: e.contact = selection.contact2 if selection.cell_um is not None: config.CELL_UM_OVERRIDE = selection.cell_um + config.ADAPTIVE_CELLS = selection.adaptive try: problem = board_io.build_problem( diff --git a/fill_resistance/quadtree.py b/fill_resistance/quadtree.py index 1126ec5..7bba038 100644 --- a/fill_resistance/quadtree.py +++ b/fill_resistance/quadtree.py @@ -123,7 +123,7 @@ def _emit(S: np.ndarray, mask: np.ndarray, max_block: int) -> LeafGrid: def _split_unbalanced(grid: LeafGrid, S: np.ndarray) -> bool: """Cap S over any leaf more than 2x an edge-adjacent neighbor, so the next emission splits it. Returns True if anything was capped.""" - ia, ib, _ = leaf_faces(grid) + ia, ib, _, _ = leaf_faces(grid) sa, sb = grid.size[ia], grid.size[ib] big = np.unique(np.concatenate([ia[sa > 2 * sb], ib[sb > 2 * sa]])) for lid in big: @@ -132,32 +132,36 @@ def _split_unbalanced(grid: LeafGrid, S: np.ndarray) -> bool: return len(big) > 0 -def leaf_faces(grid: LeafGrid) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """(ia, ib, w) for every edge-adjacent leaf pair, each pair once; - w = shared face length in fine-cell units.""" +def leaf_faces(grid: LeafGrid) -> tuple[np.ndarray, np.ndarray, np.ndarray, + np.ndarray]: + """(ia, ib, w, axis) for every edge-adjacent leaf pair, each pair + once; w = shared face length in fine-cell units. axis 0 = x-faces + (a left of b), axis 1 = y-faces (a above b).""" n = grid.n if n == 0: z = np.zeros(0, dtype=np.int64) - return z, z, z + return z, z, z, z out = [] - for sl_a, sl_b in ((np.s_[:, :-1], np.s_[:, 1:]), - (np.s_[:-1, :], np.s_[1:, :])): + for axis, (sl_a, sl_b) in enumerate(((np.s_[:, :-1], np.s_[:, 1:]), + (np.s_[:-1, :], np.s_[1:, :]))): a = grid.id_grid[sl_a].ravel().astype(np.int64) b = grid.id_grid[sl_b].ravel().astype(np.int64) ok = (a >= 0) & (b >= 0) & (a != b) key, counts = np.unique(a[ok] * n + b[ok], return_counts=True) - out.append((key // n, key % n, counts)) + out.append((key // n, key % n, counts, + np.full(len(key), axis, dtype=np.int8))) ia = np.concatenate([o[0] for o in out]) ib = np.concatenate([o[1] for o in out]) w = np.concatenate([o[2] for o in out]) - return ia, ib, w + ax = np.concatenate([o[3] for o in out]) + return ia, ib, w, ax def leaf_edges(grid: LeafGrid, sigma) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Face conductances [S]: sigma is a scalar or per-leaf array of sheet conductance. Uniform limit -> exactly sigma per face.""" - ia, ib, w = leaf_faces(grid) + ia, ib, w, _ = leaf_faces(grid) sig = np.broadcast_to(np.asarray(sigma, dtype=float), (grid.n,)) g = w / (grid.size[ia] / (2.0 * sig[ia]) + grid.size[ib] / (2.0 * sig[ib])) @@ -166,6 +170,6 @@ def leaf_edges(grid: LeafGrid, sigma) -> tuple[np.ndarray, np.ndarray, def balanced(grid: LeafGrid) -> bool: """2:1 balance invariant: edge-adjacent leaves differ <= 2x in size.""" - ia, ib, _ = leaf_faces(grid) + ia, ib, _, _ = leaf_faces(grid) sa, sb = grid.size[ia], grid.size[ib] return bool((np.maximum(sa, sb) <= 2 * np.minimum(sa, sb)).all()) diff --git a/fill_resistance/solver.py b/fill_resistance/solver.py index 49ce5ed..4a14982 100644 --- a/fill_resistance/solver.py +++ b/fill_resistance/solver.py @@ -134,6 +134,56 @@ def _sigma_2d(stack: RasterStack, li: int, sigma_layer: float, return s +def _barrel_links(stack: RasterStack, problem: Problem + ) -> tuple[list, int]: + """Vertical barrel links on the fine grid, shared by the uniform and + adaptive paths: [(via_index, layer_a, i_a, j_a, layer_b, i_b, j_b, + r_dc), ...] plus the count of dead barrels (span >= 2 layers but + reached copper on < 2). Connection cell per layer: the cell under + the barrel, or the nearest copper cell whose center lies within the + pad footprint (+1 cell of rasterization slop) - fills joined to the + barrel by thermal-relief spokes still connect, wider antipads do not + (the barrel then bridges the layers above/below).""" + L, ny, nx = stack.masks.shape + h = stack.h_nm + links = [] + dead = 0 + for vi, via in enumerate(problem.vias): + cell = stack.cell_of(via.x, via.y) + if cell is None: + continue + i, j = cell + span = [li for li, layer in enumerate(problem.layers) + if via.spans(layer.z_nm)] + r_nm = max(via.pad_nm, via.drill_nm + 300_000) / 2.0 + h + win = int(r_nm // h) + 1 + i0, i1 = max(0, i - win), min(ny, i + win + 1) + j0, j1 = max(0, j - win), min(nx, j + win + 1) + xs = stack.x0_nm + (np.arange(j0, j1) + 0.5) * h - via.x + ys = stack.y0_nm + (np.arange(i0, i1) + 0.5) * h - via.y + d2 = ys[:, None] ** 2 + xs[None, :] ** 2 + d2 = np.where(d2 <= r_nm * r_nm, d2, np.inf) + present = [] # (layer, i, j) per layer + for li in span: + if stack.masks[li, i, j]: + present.append((li, i, j)) + continue + dc = np.where(stack.masks[li, i0:i1, j0:j1], d2, np.inf) + ci, cj = np.unravel_index(int(np.argmin(dc)), dc.shape) + if np.isfinite(dc[ci, cj]): + present.append((li, i0 + ci, j0 + cj)) + if len(span) >= 2 and len(present) < 2: + dead += 1 + for (la, ia, ja), (lb, ib, jb) in zip(present[:-1], present[1:]): + length = problem.layers[lb].z_nm - problem.layers[la].z_nm + if length <= 0: + continue + r_dc = via.barrel_resistance(length, problem.rho_ohm_m, + problem.plating_nm) + links.append((vi, la, ia, ja, lb, ib, jb, r_dc)) + return links, dead + + def build_edges(stack: RasterStack, problem: Problem, sigmas: list[float], via_factor: float = 1.0, sigma_buildup: float = 0.0) -> Edges: @@ -175,49 +225,12 @@ def build_edges(stack: RasterStack, problem: Problem, sigmas: list[float], ww.append(2.0 * s_a * s_b / (s_a + s_b)) vv.append(np.full(len(a), -1, dtype=np.int32)) - h = stack.h_nm - dead_barrels = 0 - for vi, via in enumerate(problem.vias): - cell = stack.cell_of(via.x, via.y) - if cell is None: - continue - i, j = cell - span = [li for li, layer in enumerate(problem.layers) - if via.spans(layer.z_nm)] - # Connection cell per layer: the cell under the barrel, or the - # nearest copper cell whose center lies within the pad footprint - # (+1 cell of rasterization slop) - fills joined to the barrel by - # thermal-relief spokes still connect, wider antipads do not (the - # barrel then bridges the layers above/below as before). - r_nm = max(via.pad_nm, via.drill_nm + 300_000) / 2.0 + h - win = int(r_nm // h) + 1 - i0, i1 = max(0, i - win), min(ny, i + win + 1) - j0, j1 = max(0, j - win), min(nx, j + win + 1) - xs = stack.x0_nm + (np.arange(j0, j1) + 0.5) * h - via.x - ys = stack.y0_nm + (np.arange(i0, i1) + 0.5) * h - via.y - d2 = ys[:, None] ** 2 + xs[None, :] ** 2 - d2 = np.where(d2 <= r_nm * r_nm, d2, np.inf) - present = [] # (layer, i, j) per layer - for li in span: - if stack.masks[li, i, j]: - present.append((li, i, j)) - continue - dc = np.where(stack.masks[li, i0:i1, j0:j1], d2, np.inf) - ci, cj = np.unravel_index(int(np.argmin(dc)), dc.shape) - if np.isfinite(dc[ci, cj]): - present.append((li, i0 + ci, j0 + cj)) - if len(span) >= 2 and len(present) < 2: - dead_barrels += 1 - for (la, ia, ja), (lb, ib, jb) in zip(present[:-1], present[1:]): - length = problem.layers[lb].z_nm - problem.layers[la].z_nm - if length <= 0: - continue - r = via.barrel_resistance(length, problem.rho_ohm_m, - problem.plating_nm) * via_factor - aa.append(np.array([la * plane + ia * nx + ja], dtype=np.int64)) - bb.append(np.array([lb * plane + ib * nx + jb], dtype=np.int64)) - ww.append(np.array([1.0 / r])) - vv.append(np.array([vi], dtype=np.int32)) + links, dead_barrels = _barrel_links(stack, problem) + for vi, la, ia, ja, lb, ib, jb, r_dc in links: + aa.append(np.array([la * plane + ia * nx + ja], dtype=np.int64)) + bb.append(np.array([lb * plane + ib * nx + jb], dtype=np.int64)) + ww.append(np.array([1.0 / (r_dc * via_factor)])) + vv.append(np.array([vi], dtype=np.int32)) if stack.chain_edges is not None and len(stack.chain_edges[0]): ca, cb, cg, cl = stack.chain_edges @@ -409,22 +422,11 @@ def _face_current_density(V2: np.ndarray, mask2: np.ndarray, sigma: float, return Jmag -def _solve_equipotential(stack, e1, e2, edges): - """Dirichlet terminals at 1 V / 0 V. Returns (Vflat_unit, R, I1, I2, - mismatch, volts_per_amp, info). Fields at 1 V drive; scale by - i_test * R to get volts at I_test.""" - if (layer := electrodes_touch(stack, e1, e2)) is not None: - raise ElectrodeError( - f"The terminals touch on {layer}. With the equipotential " - f"contact model at least one cell of copper must separate " - f"them; the uniform-injection model allows touching contacts." - ) - state = np.zeros(stack.masks.size, dtype=np.uint8) - state[stack.masks.ravel()] = 1 - state[e1.ravel()] = 2 - state[e2.ravel()] = 3 - - A, rhs, idx = _assemble(state, edges, None) +def _equipotential_core(state: np.ndarray, edges: Edges): + """Dirichlet solve on any node space (fine cells or leaves): state + codes 0 off / 1 free / 2 V+ / 3 V-. Returns (Vflat_unit, R, I1, I2, + mismatch, volts_per_amp, info); fields at 1 V drive.""" + A, rhs, _ = _assemble(state, edges, None) x, info = solve_system(A, rhs) Vflat = np.zeros(state.size) @@ -440,29 +442,14 @@ def _solve_equipotential(stack, e1, e2, edges): return Vflat, R, I1, I2, mismatch, R, info -def _solve_uniform(stack, e1, e2, edges): - """Uniform orthogonal injection: every contact cell sources (sinks) - 1 A / N. Grounded at one V- cell. Returns like _solve_equipotential; - fields are at 1 A drive, so volts_per_amp = 1.""" - n = stack.masks.size - e1f, e2f = e1.ravel(), e2.ravel() - n1, n2 = int(e1f.sum()), int(e2f.sum()) - - inj = np.zeros(n) - inj[e1f] = 1.0 / n1 - inj[e2f] = -1.0 / n2 - - state = np.zeros(n, dtype=np.uint8) - state[stack.masks.ravel()] = 1 - ground = int(np.flatnonzero(e2f)[0]) - state[ground] = 3 # single Dirichlet 0 V reference; - # its sink share is exactly the flux that exits through the reference, - # so the grounded solution equals the pure-Neumann one - - A, rhs, idx = _assemble(state, edges, inj) +def _uniform_core(state: np.ndarray, inj: np.ndarray, e1f: np.ndarray, + e2f: np.ndarray, edges: Edges): + """Uniform-injection solve on any node space. state must carry the + single ground node (code 3); inj the per-node current shares.""" + A, rhs, _ = _assemble(state, edges, inj) x, info = solve_system(A, rhs) - Vflat = np.zeros(n) + Vflat = np.zeros(state.size) Vflat[state == 1] = x v_plus = float(Vflat[e1f].mean()) @@ -478,6 +465,42 @@ def _solve_uniform(stack, e1, e2, edges): return Vflat, R, 1.0, 1.0, res, 1.0, info +def _solve_equipotential(stack, e1, e2, edges): + """Dirichlet terminals at 1 V / 0 V on the uniform grid. Fields at + 1 V drive; scale by i_test * R to get volts at I_test.""" + if (layer := electrodes_touch(stack, e1, e2)) is not None: + raise ElectrodeError( + f"The terminals touch on {layer}. With the equipotential " + f"contact model at least one cell of copper must separate " + f"them; the uniform-injection model allows touching contacts." + ) + state = np.zeros(stack.masks.size, dtype=np.uint8) + state[stack.masks.ravel()] = 1 + state[e1.ravel()] = 2 + state[e2.ravel()] = 3 + return _equipotential_core(state, edges) + + +def _solve_uniform(stack, e1, e2, edges): + """Uniform orthogonal injection on the uniform grid: every contact + cell sources (sinks) 1 A / N, grounded at one V- cell (its sink + share is exactly the flux that exits through the reference, so the + grounded solution equals the pure-Neumann one). Fields at 1 A.""" + n = stack.masks.size + e1f, e2f = e1.ravel(), e2.ravel() + n1, n2 = int(e1f.sum()), int(e2f.sum()) + + inj = np.zeros(n) + inj[e1f] = 1.0 / n1 + inj[e2f] = -1.0 / n2 + + state = np.zeros(n, dtype=np.uint8) + state[stack.masks.ravel()] = 1 + ground = int(np.flatnonzero(e2f)[0]) + state[ground] = 3 + return _uniform_core(state, inj, e1f, e2f, edges) + + def _part_currents(parts, Ie, edges, e_flat, scale, i_test, contact_model, n_terminal_cells): """Current through each contact part @ I_test. Equipotential: exact @@ -498,18 +521,12 @@ def _part_currents(parts, Ie, edges, e_flat, scale, return out -def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray, - e2: np.ndarray, i_test: float, freq_hz: float = 0.0, - contact_model: str | None = None, - parts1: list | None = None, - parts2: list | None = None) -> Result: - timings = {} - L, ny, nx = stack.masks.shape - h_m = stack.h_nm * 1e-9 - if contact_model is None: - contact_model = config.CONTACT_MODEL - - # effective (AC) sheet conductances and barrel factor +def _conductance_params(problem: Problem, stack: RasterStack, + freq_hz: float): + """Effective (possibly AC) sheet conductances per layer, Rs ratios, + barrel factor and buildup conductance - shared by the uniform-grid + and adaptive solve paths.""" + L = stack.nlayers sigmas = [ 1.0 / skin.sheet_resistance_ac( problem.layers[li].thickness_nm * 1e-9, freq_hz, @@ -544,6 +561,28 @@ def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray, f"per-layer Rs ratio " f"{', '.join(f'{r:.2f}' for r in rs_ratios)}, " f"via factor {via_factor:.2f}") + return sigmas, rs_ratios, via_factor, sigma_buildup + + +def run_solve(problem: Problem, stack: RasterStack, e1: np.ndarray, + e2: np.ndarray, i_test: float, freq_hz: float = 0.0, + contact_model: str | None = None, + parts1: list | None = None, + parts2: list | None = None) -> Result: + if contact_model is None: + contact_model = config.CONTACT_MODEL + if config.ADAPTIVE_CELLS: + from . import adaptive + return adaptive.run_solve_adaptive(problem, stack, e1, e2, i_test, + freq_hz, contact_model, + parts1, parts2) + + timings = {} + L, ny, nx = stack.masks.shape + h_m = stack.h_nm * 1e-9 + + 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) diff --git a/fill_resistance/standalone.py b/fill_resistance/standalone.py index f27c13d..d43569f 100644 --- a/fill_resistance/standalone.py +++ b/fill_resistance/standalone.py @@ -47,6 +47,9 @@ def main(argv=None) -> int: ap.add_argument("--force-iterative", action="store_true", help="use the iterative solver (AMG-CG, or Jacobi-CG " "without pyamg) regardless of problem size") + ap.add_argument("--adaptive", action="store_true", + help="solve on the adaptive quadtree grid (coarse " + "plane interiors)") args = ap.parse_args(argv) if args.cell_um is not None: @@ -55,6 +58,8 @@ def main(argv=None) -> int: config.INTERACTIVE = False if args.force_iterative: config.SPSOLVE_MAX_UNKNOWNS = 0 + if args.adaptive: + config.ADAPTIVE_CELLS = True problem = load_problem(args.dump) if args.strip_buildup: diff --git a/tests/test_adaptive.py b/tests/test_adaptive.py new file mode 100644 index 0000000..3b12cae --- /dev/null +++ b/tests/test_adaptive.py @@ -0,0 +1,156 @@ +"""Adaptive solve path (phase 2): full run_solve equivalence against the +uniform grid across the feature set. On piecewise-linear fields (strips) +the leaf system is EXACT, so those compare at solver precision.""" +import numpy as np +import pytest + +from fill_resistance import config, raster, solver +from fill_resistance.geometry import Electrode, LayerFill, Polygon, Problem, \ + TrackSeg +from tests.test_capping import _two_layer +from tests.util import NM, make_multilayer, make_problem, rect_mm, \ + strip_problem + + +def _run(problem, h_mm, model="equipotential", adaptive=False, + monkeypatch=None, parts=False, freq=0.0): + if monkeypatch is not None: + monkeypatch.setattr(config, "ADAPTIVE_CELLS", adaptive) + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + kw = {} + if parts: + p1, p2 = raster.electrode_partition(stack, problem) + kw = dict(parts1=p1, parts2=p2) + return solver.run_solve(problem, stack, e1, e2, 1.0, freq, + contact_model=model, **kw) + + +def test_strip_close_both_models(monkeypatch): + """Uniform strip, both contact models. Coarse-fine interfaces carry + a first-order tangential flux error (laterally offset leaf centers), + so adaptive R sits up to ~2% LOW of the production R - the narrow + strip is the worst case (transition rings span most of the width).""" + for model in ("equipotential", "uniform"): + p = strip_problem(length=50, width=10, e_len=5) + ref = _run(p, 0.25, model, adaptive=False, monkeypatch=monkeypatch) + p2 = strip_problem(length=50, width=10, e_len=5) + ada = _run(p2, 0.25, model, adaptive=True, monkeypatch=monkeypatch) + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=0.02), model + assert ada.R_ohm <= ref.R_ohm * 1.001 # bias is low, not high + assert ada.n_free < ref.n_free + assert ada.power_balance_rel < 1e-9 + + +def test_plate_with_holes_close(monkeypatch): + holes = [] + for i in range(5): + for j in range(5): + x, y = 8 * i + 3, 8 * j + 3 + holes.append([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)]) + outline = [(0, 0), (40, 0), (40, 40), (0, 40)] + + def prob(): + return make_problem([(outline, holes)], rect1_mm=(0, 15, 2, 25), + rect2_mm=(38, 15, 40, 25)) + + ref = _run(prob(), 0.1, adaptive=False, monkeypatch=monkeypatch) + ada = _run(prob(), 0.1, adaptive=True, monkeypatch=monkeypatch) + assert ada.n_free < 0.5 * ref.n_free + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=0.01) + + +def test_via_chain_exact(monkeypatch): + """1-cell strips + via: everything is keep-fine or boundary, so the + adaptive path must reproduce the exact discrete solution.""" + STRIP = [(0, 0), (10, 0), (10, 1), (0, 1)] + + def prob(): + return make_multilayer( + [[(STRIP, [])], [(STRIP, [])]], + rect1_mm=(0, 0, 1, 1), rect2_mm=(9, 0, 10, 1), + contact1="L0", contact2="L1", + vias_mm=[(5.5, 0.5)], gap_mm=1.0) + + ref = _run(prob(), 1.0, adaptive=False, monkeypatch=monkeypatch) + ada = _run(prob(), 1.0, adaptive=True, monkeypatch=monkeypatch) + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=1e-9) + assert len(ada.via_reports) == 1 + assert ada.via_reports[0].current_a == pytest.approx(1.0, rel=1e-9) + + +def test_1d_trace_bridge(monkeypatch): + """Sub-resolution trace bridging two pours: chain cells are pinned + fine, pours coarsen; R matches the uniform grid closely.""" + pour1 = [(0, 0), (30, 0), (30, 30), (0, 30)] + pour2 = [(50, 0), (80, 0), (80, 30), (50, 30)] + + def prob(): + seg = TrackSeg(layer_name="F.Cu", + points=np.array([[15 * NM, 15 * NM], + [65 * NM, 15 * NM]], dtype=np.int64), + width_nm=int(0.2 * NM)) + return Problem( + board_path="synthetic", net_name="TEST", rho_ohm_m=1.68e-8, + plating_nm=18_000, + layers=[LayerFill( + layer_name="F.Cu", thickness_nm=70_000, z_nm=0, + polygons=[Polygon(outline=(np.array(pour1) * NM + ).astype(np.int64)), + Polygon(outline=(np.array(pour2) * NM + ).astype(np.int64))])], + vias=[], + electrodes1=[Electrode(rect=rect_mm((0, 10, 2, 20)))], + electrodes2=[Electrode(rect=rect_mm((78, 10, 80, 20)))], + tracks=[seg]) + + ref = _run(prob(), 0.5, adaptive=False, monkeypatch=monkeypatch) + ada = _run(prob(), 0.5, adaptive=True, monkeypatch=monkeypatch) + assert ada.n_free < 0.6 * ref.n_free + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=2e-3) + + +def test_buildup_close_on_strip(monkeypatch): + """Half-coverage buildup strip: buildup cells are pinned fine; the + remaining deviation is the interface flux bias of the plain half.""" + from tests.test_buildup import _with_buildup + + def prob(): + return _with_buildup(strip_problem(length=50, width=10, e_len=5), + [[(25, 0), (50, 0), (50, 10), (25, 10)]]) + + ref = _run(prob(), 0.5, adaptive=False, monkeypatch=monkeypatch) + ada = _run(prob(), 0.5, adaptive=True, monkeypatch=monkeypatch) + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=0.02) + + +def test_capped_via_close(monkeypatch): + """Ring + thin-cap mouth (thick_scale) under the adaptive grid.""" + ref = _run(_two_layer(drill_mm=2.0, pad_mm=2.6), 0.25, + adaptive=False, monkeypatch=monkeypatch) + ada = _run(_two_layer(drill_mm=2.0, pad_mm=2.6), 0.25, + adaptive=True, monkeypatch=monkeypatch) + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=2e-3) + + +def test_part_currents_and_ac(monkeypatch): + """Per-part currents and the AC path work on leaves.""" + p = strip_problem(length=50, width=10, e_len=5) + ref = _run(p, 0.5, adaptive=False, monkeypatch=monkeypatch, + parts=True, freq=2e6) + p2 = strip_problem(length=50, width=10, e_len=5) + ada = _run(p2, 0.5, adaptive=True, monkeypatch=monkeypatch, + parts=True, freq=2e6) + assert ada.R_ohm == pytest.approx(ref.R_ohm, rel=0.02) + assert ada.part_currents1[0][1] == pytest.approx( + ref.part_currents1[0][1], rel=1e-9) # single part = full current + assert ada.rs_ratios == ref.rs_ratios + + +def test_max_cell_size_respected(monkeypatch): + from fill_resistance.adaptive import _max_block + assert _max_block(100_000.0) == 16 # 2000 um / 100 um cells + monkeypatch.setattr(config, "ADAPTIVE_MAX_CELL_UM", 250.0) + assert _max_block(100_000.0) == 2 + monkeypatch.setattr(config, "ADAPTIVE_MAX_CELL_UM", 50.0) + assert _max_block(100_000.0) == 1 # never below the fine cell \ No newline at end of file