commit 06c62e04f8593b0b4228cc87b72211c26b1eee5c Author: janik Date: Tue Jul 14 17:22:00 2026 +0700 Initial import: KiCad zone resistance plugin DC/AC resistance, power dissipation, and via/injection-area currents of copper zone fills. KiCad 10 IPC-API plugin (kicad-python/kipy): multi-layer via-coupled FDM solver, multi-part terminals via User.1/User.2 marker layers, pads as contacts, uniform-injection and equipotential contact models, per-foil skin effect, optional solder/copper buildup on mask openings. 54-case test suite incl. exact analytic references. Co-Authored-By: Claude Fable 5 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..17cd2fc --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.venv/ +__pycache__/ +*.pyc +.pytest_cache/ diff --git a/README.md b/README.md new file mode 100644 index 0000000..845ae1a --- /dev/null +++ b/README.md @@ -0,0 +1,134 @@ +# Fill Resistance — KiCad 10 plugin + +Computes the **DC resistance of copper zone fills** between two contacts, +**single- or multi-layer**: the chosen net's fills on the selected copper +layers are solved as coupled finite-difference sheets linked by the net's +**via and through-hole-pad barrels** (18 µm plating, configurable). Shows +per-layer rasterized maps, potential, current density, and **power +density**, reports **per-via currents** (via ampacity!) and total +dissipation at a **selectable test current**. PNGs + a text summary are +saved per run. + +Uses the KiCad **IPC API** (`kicad-python` / `kipy`), not the deprecated +SWIG API. Requires KiCad **10.0.1+**. + +## Setup (one-time) + +1. **Enable the API server**: KiCad → Preferences → Plugins → check + *Enable KiCad API*. +2. **Check the interpreter path** on the same page: should be + `C:\Program Files\KiCad\10.0\bin\pythonw.exe` (after a 9→10 upgrade it + can point at KiCad 9). +3. **Deploy**: + ```powershell + powershell -ExecutionPolicy Bypass -File deploy.ps1 # junction (dev) + powershell -ExecutionPolicy Bypass -File deploy.ps1 -Mode Copy + ``` +4. **Restart KiCad**; first load builds the plugin venv (numpy, scipy, + matplotlib, PySide6 — takes minutes; the Ω button appears when done). + If stuck: Preferences → Plugins → *Recreate Plugin Environment*. + +## Usage + +1. Mark the current-injection terminals. Each terminal may have + **multiple parts** (all merged into one externally-bonded contact): + - **V+ rectangles on `User.1`**, **V− rectangles on `User.2`** + (marker layers, configurable via `ELECTRODE_POS_LAYER` / + `ELECTRODE_NEG_LAYER`), any number per side, axis-aligned; + - **pads** (real copper shape; through-hole pad contacts all layers, + SMD pad its own layer) — selected pads fill a side that has no + rectangles; + - legacy: exactly 2 selected contacts with no marker rectangles still + works; empty selection scans the whole board's marker layers. +2. **Select the contacts**, click the **Fill Resistance** Ω button. +3. In the **dialog**, pick the net (defaults to the selected pad's net), + check the **layers** to include, set each contact's layer scope + ("All selected layers" = bolted-lug/through contact), the **test + current**, and optionally a grid cell size. Multiple layers are coupled + through the net's via/pad barrels automatically. +4. Read R / voltage drop / total power in the figure titles and status + bar. Outputs land in `\fill_res_results\\`: + per-layer `1_raster_map` / `2_potential` / `3_current_density` / + `4_power_density` PNGs, `summary.txt` (incl. the busiest vias with + per-via current and dissipation, and the **current through each + injection area** — computed flux with the equipotential model, + prescribed area share with the uniform model), `geometry_dump.json`. + +## Model & limits + +- Sheet model per layer: R□ = ρ/t, ρ = 1.68e-8 Ωm (20 °C), t from the + board's physical stackup. Layer z-positions from the stackup drive the + barrel lengths. +- Via/pad barrels: thin-wall annulus, R = ρ·L/(π·d·t_plating), + `VIA_PLATING_UM = 18` in `fill_resistance/config.py`. Vias are always + plated; capped vs uncapped does not change the layer-to-layer DC path + (the ≥5 µm cap sits over the hole mouth in parallel with the + annular-ring contact, not in series). A barrel passing an antipad + bridges the layers above/below with the full barrel length. +- Tracks and pad copper (other than the selected contacts) are **not** + part of the conductor model — zone fills + barrels only. +- **Solder buildup on mask openings** (dialog checkbox, **off by + default**; `INCLUDE_MASK_BUILDUP`): zones drawn on `F.Mask`/`B.Mask` + are treated as mask openings that collect `SOLDER_THICKNESS_UM` + (50 µm) of solder on the exposed pour, plus an optional user-defined + added copper thickness (dialog field, e.g. a soldered busbar/wire). + The sheet conductance there becomes t_Cu/ρ_Cu + t_solder/ρ_solder + + t_extra/ρ_Cu (SAC305 ρ = 1.32e-7 Ωm: 50 µm solder ≈ 6.4 µm copper); + interface faces use harmonic-mean conductances. Buildup areas render + tin-gray on the raster map; |J| in them is referenced to the + conductance-equivalent copper thickness. +- **Contact models** (dialog / `CONTACT_MODEL`): default **uniform + injection** — a conductor pressed on top feeds the current orthogonally + with uniform surface density, so |J| ramps across the contact area + (R = ΔV̄/I from area-averaged terminal potentials); or + **equipotential** — ideal bonded lug (Dirichlet). The two bracket a + real contact: R_equipotential ≤ R_real ≤ R_uniform. +- Fields are reported at the dialog's test current; power scales with I². +- **Skin effect (f > 0)**: per-layer effective sheet resistance from the + exact 1D foil-diffusion solution `Zs = τρ·coth(τt)`, `τ = (1+j)/δ` + (`SKIN_SIDES = 1` in config: plane facing a return plane; `2` = + isolated foil), and the analogous correction for the 18 µm barrel wall. + Enter one frequency per run (e.g. a switching harmonic, with its RMS + amplitude as the test current) — suffixes `k`/`M` accepted. + **Caveat:** only through-thickness crowding is modeled. Lateral + (proximity-effect) redistribution needs a magneto-quasistatic solver + and is not captured — since the resistance-driven distribution is the + minimum-dissipation one, AC results are a rigorous **lower bound**. + Rule of thumb for 70 µm foil: skin is negligible below ~300 kHz + (δ = 173 µm at 142 kHz), ~+11 % at 1 MHz. +- 5-point FDM per layer on an auto-sized shared grid (~500 k cells total + across layers by default). Direct sparse solve up to 700 k unknowns, + Jacobi-CG above. Discretization error typically ≲ 2 % at defaults — + halve the cell size and compare to judge convergence. + +## Offline / development + +Every run writes `geometry_dump.json`; re-solve without KiCad: + +```powershell +.venv\Scripts\python.exe -m fill_resistance.standalone dump.json ` + [--current 40] [--cell-um 50] [--layers F.Cu,In1.Cu] [--no-show] ` + [--out DIR] [--force-iterative] +``` + +Dev environment, tests, headless extraction: + +```powershell +uv venv --python 3.11 .venv +uv pip install --python .venv\Scripts\python.exe kicad-python numpy scipy matplotlib pytest +.venv\Scripts\python.exe -m pytest tests -q # incl. exact analytic cases +.venv\Scripts\python.exe smoke\smoke_probe.py # IPC API probe vs live KiCad +.venv\Scripts\python.exe -m fill_resistance.board_io dump.json [NET] # extract only +``` + +## Troubleshooting + +- **No toolbar button**: venv still building (wait), or build failed → + *Recreate Plugin Environment*; check the interpreter path (setup 2). +- **"Could not connect to KiCad's IPC API"**: API server not enabled, or + KiCad not running (no headless mode in KiCad 10). +- **"KiCad is busy"**: a modal dialog is open in KiCad — close it, rerun. +- **Windows don't appear**: they may open behind KiCad (raised + best-effort); PNGs are always saved regardless. +- **Result seems too low/high**: remember the model is fills + barrels + only, with ideal contacts; measure electrode-to-electrode. diff --git a/conftest.py b/conftest.py new file mode 100644 index 0000000..2ca7a99 --- /dev/null +++ b/conftest.py @@ -0,0 +1,4 @@ +import os +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) diff --git a/deploy.ps1 b/deploy.ps1 new file mode 100644 index 0000000..53c2991 --- /dev/null +++ b/deploy.ps1 @@ -0,0 +1,39 @@ +# Deploy the plugin into the KiCad user plugins directory. +# .\deploy.ps1 -> NTFS junction (best for development) +# .\deploy.ps1 -Mode Copy -> plain copy (no repo link) +param( + [ValidateSet('Junction', 'Copy')] + [string]$Mode = 'Junction', + [string]$KiCadVersion = '10.0' +) + +$src = $PSScriptRoot +$pluginsDir = Join-Path ([Environment]::GetFolderPath('MyDocuments')) "KiCad\$KiCadVersion\plugins" +$dst = Join-Path $pluginsDir 'fill-resistance' + +if (-not (Test-Path $pluginsDir)) { + Write-Error "KiCad plugins dir not found: $pluginsDir" + exit 1 +} + +if (Test-Path $dst) { + $item = Get-Item $dst -Force + if ($item.LinkType) { + # junction: rmdir removes the link only, never the target contents + cmd /c rmdir "$dst" + } else { + Remove-Item $dst -Recurse -Force + } +} + +if ($Mode -eq 'Junction') { + New-Item -ItemType Junction -Path $dst -Target $src | Out-Null + Write-Host "junction created: $dst -> $src" +} else { + $exclude = @('.venv', '.git', 'tests', 'smoke', '__pycache__', '.pytest_cache') + New-Item -ItemType Directory -Force $dst | Out-Null + Get-ChildItem $src -Force | Where-Object { $exclude -notcontains $_.Name } | + ForEach-Object { Copy-Item $_.FullName -Destination $dst -Recurse -Force } + Write-Host "copied plugin to: $dst" +} +Write-Host "Restart KiCad (or refresh plugins) and wait for the plugin venv build." diff --git a/fill_res_action.py b/fill_res_action.py new file mode 100644 index 0000000..d040a73 --- /dev/null +++ b/fill_res_action.py @@ -0,0 +1,14 @@ +"""Entrypoint launched by KiCad when the toolbar button is pressed. + +Lives at the plugin root (next to plugin.json) and only bootstraps the +package import path, so fill_resistance/ can use normal absolute imports. +""" +import os +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) + +from fill_resistance.main import main + +if __name__ == "__main__": + main() diff --git a/fill_resistance/__init__.py b/fill_resistance/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/fill_resistance/board_io.py b/fill_resistance/board_io.py new file mode 100644 index 0000000..fa67512 --- /dev/null +++ b/fill_resistance/board_io.py @@ -0,0 +1,480 @@ +"""All KiCad IPC access. This is the ONLY module that imports kipy; +everything downstream works on plain geometry dataclasses. + +Run `python -m fill_resistance.board_io dump.json [net]` against a live +KiCad to extract without the dialog (all layers of the net, defaults). +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path + +from kipy import KiCad +from kipy.board import Board +from kipy.board_types import BoardRectangle, Pad +from kipy.proto.board.board_pb2 import BoardStackupLayerType +from kipy.proto.board.board_types_pb2 import ZoneType +from kipy.util.board_layer import (canonical_name, is_copper_layer, + layer_from_canonical_name) + +from . import config +from .errors import ApiVersionError, CandidateError, SelectionError +from .geometry import (Electrode, LayerFill, Polygon, Problem, Rect, + SurfaceBuildup, ViaLink, linearize_ring) + +MASK_TO_COPPER = {"F.Mask": "F.Cu", "B.Mask": "B.Cu"} + +# zone fills are polygonal in practice; tolerance only guards arc nodes +ARC_TOL_NM = 10_000 + + +def connect() -> tuple[KiCad, Board]: + try: + kicad = KiCad() + kicad.ping() + except Exception as e: + raise ApiVersionError( + f"Could not connect to KiCad's IPC API: {e}\n" + f"Is KiCad running with the API server enabled " + f"(Preferences > Plugins > Enable KiCad API)?" + ) + try: + print(f"connected to KiCad {kicad.get_version()}") + except Exception: + pass + try: + board = kicad.get_board() + except Exception as e: + raise SelectionError( + f"Could not get the open board from KiCad: {e}\n" + f"Open the PCB in the board editor and run again." + ) + return kicad, board + + +def board_dir(board: Board) -> Path: + # document.board_filename is a bare file name (no directory) in + # KiCad 10.0.1; the project path is the reliable location + try: + path = board.get_project().path + if path and Path(path).is_dir(): + return Path(path) + except Exception: + pass + try: + filename = getattr(board.document, "board_filename", "") or "" + if Path(filename).is_absolute(): + return Path(filename).parent + except Exception: + pass + return Path.cwd() + + +# --- stackup geometry -------------------------------------------------------- + +@dataclass +class StackupInfo: + names: list[str] # copper layers, top to bottom + thickness_nm: dict[str, int] + z_nm: dict[str, int] # copper center depth + z_bot_nm: int # total stack thickness + + +def get_stackup_info(board: Board) -> StackupInfo: + names: list[str] = [] + thickness: dict[str, int] = {} + z_center: dict[str, int] = {} + z = 0 + for sl in board.get_stackup().layers: + t = int(sl.thickness or 0) + if sl.type == BoardStackupLayerType.BSLT_COPPER: + name = canonical_name(sl.layer) + if t <= 0: + t = int(config.FALLBACK_THICKNESS_UM * 1000) + print(f"warning: stackup gives no thickness for {name}; " + f"assuming {config.FALLBACK_THICKNESS_UM} um") + names.append(name) + thickness[name] = t + z_center[name] = z + t // 2 + z += t + if not names: + raise CandidateError( + "Could not read any copper layer from the board stackup." + ) + return StackupInfo(names=names, thickness_nm=thickness, z_nm=z_center, + z_bot_nm=z) + + +# --- electrodes from selection ---------------------------------------------- + +def _box2_to_rect(box, layer_name: str) -> Rect: + try: + pos, size = box.pos, box.size + return Rect.normalized(pos.x, pos.y, pos.x + size.x, pos.y + size.y, + layer_name) + except AttributeError: + c, s = box.center, box.size + return Rect.normalized(c.x - s.x // 2, c.y - s.y // 2, + c.x + s.x // 2, c.y + s.y // 2, layer_name) + + +def _convert_poly(poly_with_holes) -> Polygon: + def ring(polyline): + nodes = [] + for node in polyline.nodes: + if node.has_point: + nodes.append(("pt", (node.point.x, node.point.y))) + elif node.has_arc: + arc = node.arc + nodes.append(("arc", ((arc.start.x, arc.start.y), + (arc.mid.x, arc.mid.y), + (arc.end.x, arc.end.y)))) + return linearize_ring(nodes, ARC_TOL_NM) + + return Polygon(outline=ring(poly_with_holes.outline), + holes=[ring(h) for h in poly_with_holes.holes]) + + +def _pad_drill_nm(pad_or_via) -> int: + try: + return int(pad_or_via.padstack.drill.diameter.x) + except Exception: + return 0 + + +def _pad_default_contact(pad: Pad) -> str: + if _pad_drill_nm(pad) > 0: + return "all" # through-hole: contacts the stack + try: + copper = [canonical_name(l) for l in pad.padstack.layers + if is_copper_layer(l)] + if len(copper) == 1: + return copper[0] # SMD: its own layer + except Exception: + pass + return "all" + + +def _pad_polygons(board: Board, pad: Pad, contact: str) -> list[Polygon] | None: + layer_ids = [] + if contact != "all": + try: + layer_ids.append(layer_from_canonical_name(contact)) + except Exception: + pass + for name in ("F.Cu", "B.Cu"): + try: + layer_ids.append(layer_from_canonical_name(name)) + except Exception: + pass + for lid in layer_ids: + try: + shape = board.get_pad_shapes_as_polygons(pad, layer=lid) + if shape is not None: + return [_convert_poly(shape)] + except Exception: + continue + return None + + +def _to_electrode(board: Board, item) -> Electrode: + if isinstance(item, BoardRectangle): + tl, br = item.top_left, item.bottom_right + rect = Rect.normalized(tl.x, tl.y, br.x, br.y, + canonical_name(item.layer)) + cx = (rect.x0 + rect.x1) / 2e6 + cy = (rect.y0 + rect.y1) / 2e6 + return Electrode(rect=rect, contact="all", + label=f"rect({cx:.1f},{cy:.1f})") + # Pad + pad: Pad = item + contact = _pad_default_contact(pad) + net = pad.net.name if pad.net is not None else "?" + label = f"pad {pad.number}@{net}" + box = board.get_item_bounding_box(pad) + if box is None: + raise SelectionError(f"Could not get the bounding box of {label}.") + rect = _box2_to_rect(box, "pad") + return Electrode(rect=rect, contact=contact, + polygons=_pad_polygons(board, pad, contact), label=label) + + +def _net_hint_of(pads: list[Pad]) -> str | None: + for pad in pads: + if pad.net is not None: + return pad.net.name + return None + + +def get_electrodes(board: Board + ) -> tuple[list[Electrode], list[Electrode], str | None]: + """Terminals from the selection. Each terminal may have MULTIPLE parts + (all merged into one externally-bonded contact): + + - rectangles on ELECTRODE_POS_LAYER -> V+ parts, on ELECTRODE_NEG_LAYER + -> V- parts; selected pads fill a side that has no rectangles; + - no marker rectangles selected: legacy mode, exactly 2 items + (rects/pads, any layer) -> one part each; + - empty selection: board-wide scan of both marker layers. + """ + pos_l = config.ELECTRODE_POS_LAYER + neg_l = config.ELECTRODE_NEG_LAYER + scheme = (f"Draw V+ rectangle(s) on {pos_l} and V- rectangle(s) on " + f"{neg_l} (axis-aligned), and/or select pads for a side " + f"without rectangles.") + + selection = list(board.get_selection()) + rects = [s for s in selection if isinstance(s, BoardRectangle)] + pads = [s for s in selection if isinstance(s, Pad)] + + if not selection: + allr = [s for s in board.get_shapes() if isinstance(s, BoardRectangle)] + pos = [r for r in allr if canonical_name(r.layer) == pos_l] + neg = [r for r in allr if canonical_name(r.layer) == neg_l] + if pos and neg: + print(f"selection empty - using {len(pos)} rectangle(s) on " + f"{pos_l} as V+ and {len(neg)} on {neg_l} as V-") + return ([_to_electrode(board, r) for r in pos], + [_to_electrode(board, r) for r in neg], None) + raise SelectionError( + f"Nothing selected, and the board-wide scan found " + f"{len(pos)} rectangle(s) on {pos_l} / {len(neg)} on {neg_l} " + f"(need at least one on each).\n{scheme}" + ) + + pos = [r for r in rects if canonical_name(r.layer) == pos_l] + neg = [r for r in rects if canonical_name(r.layer) == neg_l] + other = [r for r in rects if canonical_name(r.layer) not in (pos_l, neg_l)] + + if pos or neg: + if other: + raise SelectionError( + f"{len(other)} selected rectangle(s) are on neither marker " + f"layer ({pos_l} = V+, {neg_l} = V-). {scheme}" + ) + es1 = [_to_electrode(board, r) for r in pos] + es2 = [_to_electrode(board, r) for r in neg] + if pads and es1 and es2: + raise SelectionError( + f"Cannot assign the {len(pads)} selected pad(s): both marker " + f"layers already provide rectangles. Use pads only for a " + f"side that has none." + ) + if pads: + pad_parts = [_to_electrode(board, p) for p in pads] + if not es1: + es1 = pad_parts + else: + es2 = pad_parts + if es1 and es2: + return es1, es2, _net_hint_of(pads) + raise SelectionError( + f"Only one terminal defined: V+ has {len(es1)} and V- has " + f"{len(es2)} contact(s). {scheme}" + ) + + items = rects + pads + if len(items) == 2: + return ([_to_electrode(board, items[0])], + [_to_electrode(board, items[1])], _net_hint_of(pads)) + raise SelectionError( + f"The selection has {len(rects)} rectangle(s) (none on the marker " + f"layers) and {len(pads)} pad(s); without marker layers exactly 2 " + f"contacts are needed.\n{scheme}" + ) + + +# --- fills ------------------------------------------------------------------- + +def gather_net_fills(board: Board) -> dict[str, dict[str, list[Polygon]]]: + """net -> layer_name -> merged fill polygons (non-empty only).""" + fills: dict[str, dict[str, list[Polygon]]] = {} + for zone in board.get_zones(): + if zone.type != ZoneType.ZT_COPPER: + continue + net = zone.net.name if zone.net is not None else "" + for layer, polys in zone.filled_polygons.items(): + if not is_copper_layer(layer) or not polys: + continue + fills.setdefault(net, {}).setdefault( + canonical_name(layer), []).extend( + _convert_poly(p) for p in polys) + return fills + + +def _rect_overlaps(rect: Rect, polygons: list[Polygon]) -> bool: + for p in polygons: + px0, py0 = p.outline.min(axis=0) + px1, py1 = p.outline.max(axis=0) + if rect.x0 <= px1 and rect.x1 >= px0 and rect.y0 <= py1 and rect.y1 >= py0: + return True + return False + + +def nets_overlapping(fills: dict, es1: list[Electrode], + es2: list[Electrode]) -> list[str]: + """Nets whose fills overlap both terminals (any part, any layer each - + the connection may go through vias). Permissive bbox prefilter.""" + out = [] + for net, per_layer in fills.items(): + hit1 = any(_rect_overlaps(e.rect, polys) for e in es1 + for polys in per_layer.values()) + hit2 = any(_rect_overlaps(e.rect, polys) for e in es2 + for polys in per_layer.values()) + if hit1 and hit2: + out.append(net) + return sorted(out) + + +def gather_mask_buildups(board: Board) -> dict[str, list[Polygon]]: + """Zones on F.Mask/B.Mask (mask openings) -> fill polygons keyed by + the outer copper layer they expose.""" + out: dict[str, list[Polygon]] = {} + for zone in board.get_zones(): + try: + filled = zone.filled_polygons + except Exception: + continue + for layer, polys in filled.items(): + copper = MASK_TO_COPPER.get(canonical_name(layer)) + if copper and polys: + out.setdefault(copper, []).extend( + _convert_poly(p) for p in polys) + return out + + +def any_zone_unfilled(board: Board) -> bool: + return any(z.type == ZoneType.ZT_COPPER and not z.filled + for z in board.get_zones()) + + +def refill(board: Board) -> None: + print("refilling zones - this modifies the open document ...") + board.refill_zones(block=True) + + +# --- barrels ----------------------------------------------------------------- + +def _padstack_span(padstack, stackup: StackupInfo) -> tuple[int, int]: + """(z_top, z_bot) of the barrel; falls back to the full stack.""" + try: + copper = [canonical_name(l) for l in padstack.layers + if is_copper_layer(l)] + zs = [stackup.z_nm[c] for c in copper if c in stackup.z_nm] + if len(zs) >= 2: + return min(zs) - 1, max(zs) + 1 + except Exception: + pass + return -1, stackup.z_bot_nm + 1 + + +def gather_barrels(board: Board, net_name: str, + stackup: StackupInfo) -> list[ViaLink]: + barrels = [] + for via in board.get_vias(): + if via.net is None or via.net.name != net_name: + continue + drill = int(via.drill_diameter or 0) or _pad_drill_nm(via) + if drill <= 0: + continue + z_top, z_bot = _padstack_span(via.padstack, stackup) + barrels.append(ViaLink(x=via.position.x, y=via.position.y, + drill_nm=drill, z_top_nm=z_top, + z_bot_nm=z_bot, kind="via")) + if config.INCLUDE_TH_PADS: + for pad in board.get_pads(): + if pad.net is None or pad.net.name != net_name: + continue + drill = _pad_drill_nm(pad) + if drill <= 0: + continue + barrels.append(ViaLink(x=pad.position.x, y=pad.position.y, + drill_nm=drill, z_top_nm=-1, + z_bot_nm=stackup.z_bot_nm + 1, kind="pad")) + return barrels + + +# --- top level ---------------------------------------------------------------- + +def build_problem(board: Board, net: str, layer_names: list[str], + es1: list[Electrode], es2: list[Electrode], + stackup: StackupInfo, fills: dict, + buildups: dict[str, list[Polygon]] | None = None, + extra_cu_um: float | None = None) -> Problem: + per_layer = fills.get(net, {}) + layers = [] + for name in stackup.names: # keep stackup order + if name not in layer_names: + continue + polys = per_layer.get(name, []) + if not polys: + print(f"note: net {net} has no fill on {name} - layer skipped") + continue + if config.COPPER_THICKNESS_UM is not None: + t, source = int(config.COPPER_THICKNESS_UM * 1000), "override" + else: + t, source = stackup.thickness_nm[name], "stackup" + layers.append(LayerFill(layer_name=name, thickness_nm=t, + z_nm=stackup.z_nm[name], polygons=polys)) + if not layers: + raise CandidateError( + f"Net {net} has no fill on any of the selected layers " + f"({', '.join(layer_names)})." + ) + vias = gather_barrels(board, net, stackup) if len(layers) > 1 else [] + included = {l.layer_name for l in layers} + buildup_list = [ + SurfaceBuildup(layer_name=name, polygons=polys) + for name, polys in (buildups or {}).items() if name in included + ] + print(f"net {net}: {len(layers)} layer(s) " + f"({', '.join(l.layer_name for l in layers)}), " + f"{len(vias)} via/pad barrel(s)" + + (f", solder buildup on " + f"{', '.join(b.layer_name for b in buildup_list)}" + if buildup_list else "")) + return Problem( + board_path=board.name or "", + net_name=net, + rho_ohm_m=config.RHO_CU_OHM_M, + plating_nm=int(config.VIA_PLATING_UM * 1000), + layers=layers, + vias=vias, + electrodes1=es1, + electrodes2=es2, + thickness_source=("override" if config.COPPER_THICKNESS_UM is not None + else "stackup"), + buildups=buildup_list, + solder_thickness_nm=int(config.SOLDER_THICKNESS_UM * 1000), + solder_rho_ohm_m=config.SOLDER_RHO_OHM_M, + extra_cu_nm=int((extra_cu_um if extra_cu_um is not None + else config.BUILDUP_EXTRA_CU_UM) * 1000), + ) + + +if __name__ == "__main__": + import sys + + from .geometry import save_problem + + out = Path(sys.argv[1]) if len(sys.argv) > 1 else Path("geometry_dump.json") + _, board = connect() + stackup = get_stackup_info(board) + es1, es2, net_hint = get_electrodes(board) + if any_zone_unfilled(board): + refill(board) + fills = gather_net_fills(board) + nets = nets_overlapping(fills, es1, es2) + if len(sys.argv) > 2: + net = sys.argv[2] + elif net_hint in nets: + net = net_hint + elif len(nets) == 1: + net = nets[0] + else: + print(f"candidate nets: {nets}; pass one as second argument") + sys.exit(1) + problem = build_problem(board, net, list(fills.get(net, {})), es1, es2, + stackup, fills) + save_problem(problem, out) + print(f"wrote {out}") diff --git a/fill_resistance/config.py b/fill_resistance/config.py new file mode 100644 index 0000000..e52b602 --- /dev/null +++ b/fill_resistance/config.py @@ -0,0 +1,72 @@ +"""All tunable constants. v1 has no GUI dialog: edit here, re-run. + +A future version may read overrides from /fill_res_config.json. +""" + +# --- Grid sizing --- +# Benchmarked on the VOUT+ plane (147x59 mm): R changes < 0.3% from +# 150 um cells down to 50 um; 1.7M unknowns direct-solve in ~17 s +# (raster ~20 s). Accuracy is feature-limited (slots/necks narrower than +# one cell), not plane-limited - override CELL_UM_OVERRIDE for boards +# with sub-cell slots. +TARGET_CELLS = 2_000_000 # auto cell size aims for roughly this many cells +HARD_MAX_CELLS = 16_000_000 # abort above this (see GridSizeError message) +MIN_CELL_UM = 25.0 # clamp for auto cell size +MAX_CELL_UM = 500.0 +CELL_UM_OVERRIDE: float | None = None # force a cell size, bypasses auto (not HARD_MAX) +MARGIN_CELLS = 2 # empty guard cells around the copper bbox + +# --- Physics --- +RHO_CU_OHM_M = 1.68e-8 # copper resistivity at 20 degC +COPPER_THICKNESS_UM: float | None = None # None -> stackup, fallback 35.0 with warning +FALLBACK_THICKNESS_UM = 35.0 +TEST_CURRENT_A = 1.0 # default injected current (dialog/CLI-selectable) +VIA_PLATING_UM = 18.0 # barrel plating thickness (always plated). + # Capped vs uncapped vias do not change the + # layer-to-layer DC path: the >=5um cap sits + # over the hole mouth in parallel with the + # annular-ring contact, not in series. +INCLUDE_TH_PADS = True # plated through-hole pads stitch layers too +SKIN_SIDES = 1 # skin-effect field config: 1 = plane facing a + # return plane (conservative), 2 = isolated foil + +# --- Solder / mask-opening buildup --- +INCLUDE_MASK_BUILDUP = False # OFF by default; dialog-toggleable. Zones on + # F.Mask/B.Mask = mask openings that collect + # solder on the pour underneath +SOLDER_THICKNESS_UM = 50.0 # solder height over opened copper +SOLDER_RHO_OHM_M = 1.32e-7 # SAC305, ~7.9x copper +BUILDUP_EXTRA_CU_UM = 0.0 # optional user-added copper (busbar/wire + # soldered into the opening); dialog-settable + +# --- Zone / layer selection --- +LAYER_HINT: str | None = None # e.g. "F.Cu" to disambiguate candidate fills +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 + +# --- Solver --- +CONTACT_MODEL = "uniform" # "uniform": conductor pressed on top injects + # orthogonally with uniform surface density + # (J ramps across the contact); "equipotential": + # ideal bonded lug (Dirichlet). The two bracket + # a real contact: R_equi <= R_real <= R_uniform. +SPSOLVE_MAX_UNKNOWNS = 2_500_000 # above this, use CG (Jacobi) instead of + # direct solve (measured: direct is ~14x + # faster at 1.7M unknowns, ~3 GB peak) +CG_TOL = 1e-8 +CG_MAXITER = 50_000 # CG iterations are cheap; large grids need many + +# --- Geometry --- +ARC_TOL_FRACTION = 0.5 # arc sagitta tolerance as a fraction of cell size + +# --- Plots / output --- +CMAP_POTENTIAL = "viridis" +CMAP_CURRENT = "inferno" +CMAP_POWER = "magma" +POWER_DYNAMIC_RANGE = 1e4 # LogNorm span for the power map (power ~ J^2) +LOG_CURRENT_SCALE = True +CURRENT_DYNAMIC_RANGE = 1e3 # LogNorm vmin = vmax / this +DPI = 150 +INTERACTIVE = True # False -> save PNGs only, never open windows +OUTPUT_DIRNAME = "fill_res_results" diff --git a/fill_resistance/dialog.py b/fill_resistance/dialog.py new file mode 100644 index 0000000..360acf6 --- /dev/null +++ b/fill_resistance/dialog.py @@ -0,0 +1,193 @@ +"""Qt selection dialog shown on plugin launch: net, layers, per-electrode +contact, test current, optional cell size. PySide6 is already a plugin +dependency (matplotlib QtAgg backend); the QApplication created here is +reused by matplotlib afterwards. +""" +from __future__ import annotations + +from dataclasses import dataclass + +from PySide6.QtCore import Qt +from PySide6.QtWidgets import (QApplication, QCheckBox, QComboBox, QDialog, + QDialogButtonBox, QFormLayout, QLabel, + QLineEdit, QListWidget, QListWidgetItem, + QVBoxLayout) + +from . import config, skin + +ALL_LAYERS = "All selected layers" +AUTO_CONTACT = "(auto: per contact part)" +MODEL_LABELS = { + "uniform": "Uniform injection (conductor pressed on top)", + "equipotential": "Equipotential (ideal bonded lug)", +} + + +@dataclass +class Selection: + net: str + layers: list[str] + contact1: str # "auto", "all" or layer name + contact2: str + current_a: float + cell_um: float | None + freq_hz: float = 0.0 + contact_model: str = "uniform" + include_buildup: bool = False + extra_cu_um: float = 0.0 + + +class _Dialog(QDialog): + def __init__(self, candidates: dict[str, list[str]], layer_order: list[str], + default_net: str, e1_label: str, e2_label: str, + contact1: str, contact2: str, buildup_layers: list[str]): + super().__init__() + self.setWindowTitle("Fill Resistance") + self.setWindowFlag(Qt.WindowStaysOnTopHint, True) + self._candidates = candidates + self._layer_order = layer_order + + form = QFormLayout() + + self.net_box = QComboBox() + for net in sorted(candidates): + self.net_box.addItem(net) + self.net_box.setCurrentText(default_net) + form.addRow("Signal (net):", self.net_box) + + self.layer_list = QListWidget() + self.layer_list.setMaximumHeight(120) + form.addRow("Layers:", self.layer_list) + + self.contact1_box = QComboBox() + self.contact2_box = QComboBox() + form.addRow(f"V+ ({e1_label}):", self.contact1_box) + form.addRow(f"V− ({e2_label}):", self.contact2_box) + + self.model_box = QComboBox() + for key in ("uniform", "equipotential"): + self.model_box.addItem(MODEL_LABELS[key], key) + default_index = 0 if config.CONTACT_MODEL == "uniform" else 1 + self.model_box.setCurrentIndex(default_index) + form.addRow("Contact model:", self.model_box) + + self.current_edit = QLineEdit(f"{config.TEST_CURRENT_A:g}") + form.addRow("Test current [A]:", self.current_edit) + + self.freq_edit = QLineEdit("") + self.freq_edit.setPlaceholderText("0 = DC (e.g. 142k, 1.5M)") + form.addRow("Frequency [Hz]:", self.freq_edit) + + self.cell_edit = QLineEdit("") + self.cell_edit.setPlaceholderText("auto") + form.addRow("Cell size [µm]:", self.cell_edit) + + self.buildup_check = QCheckBox( + f"{config.SOLDER_THICKNESS_UM:g} µm solder on mask openings" + + (f" ({', '.join(buildup_layers)})" if buildup_layers + else " (none found)")) + self.buildup_check.setChecked(bool(buildup_layers) + and config.INCLUDE_MASK_BUILDUP) + self.buildup_check.setEnabled(bool(buildup_layers)) + form.addRow("Buildup:", self.buildup_check) + + self.extracu_edit = QLineEdit(f"{config.BUILDUP_EXTRA_CU_UM:g}") + self.extracu_edit.setEnabled(bool(buildup_layers)) + form.addRow("Extra Cu in openings [µm]:", self.extracu_edit) + + buttons = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel) + buttons.accepted.connect(self.accept) + buttons.rejected.connect(self.reject) + + lay = QVBoxLayout(self) + lay.addLayout(form) + note = QLabel("Multiple layers are coupled through the net's " + "via/through-pad barrels. At f > 0 the foil-thickness " + "skin effect is applied per layer; lateral (proximity) " + "redistribution is not modeled, so AC results are a " + "lower bound.") + note.setWordWrap(True) + note.setStyleSheet("color: gray; font-size: 10px;") + lay.addWidget(note) + lay.addWidget(buttons) + + self._desired1, self._desired2 = contact1, contact2 + self.net_box.currentTextChanged.connect(self._refresh) + self._refresh() + + def _refresh(self): + net = self.net_box.currentText() + layers = [n for n in self._layer_order + if n in self._candidates.get(net, [])] + self.layer_list.clear() + for name in layers: + item = QListWidgetItem(name) + item.setFlags(item.flags() | Qt.ItemIsUserCheckable) + item.setCheckState(Qt.Checked) + self.layer_list.addItem(item) + for box, desired in ((self.contact1_box, self._desired1), + (self.contact2_box, self._desired2)): + box.clear() + box.addItem(AUTO_CONTACT) + box.addItem(ALL_LAYERS) + box.addItems(layers) + if desired == "all": + box.setCurrentText(ALL_LAYERS) + elif desired in layers: + box.setCurrentText(desired) + + def checked_layers(self) -> list[str]: + out = [] + for i in range(self.layer_list.count()): + item = self.layer_list.item(i) + if item.checkState() == Qt.Checked: + out.append(item.text()) + return out + + def selection(self) -> Selection | None: + layers = self.checked_layers() + if not layers: + return None + try: + current = float(self.current_edit.text().replace(",", ".")) + except ValueError: + current = config.TEST_CURRENT_A + cell_text = self.cell_edit.text().strip() + try: + cell = float(cell_text.replace(",", ".")) if cell_text else None + except ValueError: + cell = None + + def contact(box: QComboBox) -> str: + t = box.currentText() + if t == AUTO_CONTACT: + return "auto" + return "all" if t == ALL_LAYERS else t + + try: + extra_cu = float(self.extracu_edit.text().replace(",", ".")) + except ValueError: + extra_cu = 0.0 + return Selection(net=self.net_box.currentText(), layers=layers, + contact1=contact(self.contact1_box), + contact2=contact(self.contact2_box), + current_a=current, cell_um=cell, + freq_hz=skin.parse_frequency(self.freq_edit.text()), + contact_model=self.model_box.currentData(), + include_buildup=self.buildup_check.isChecked(), + extra_cu_um=max(0.0, extra_cu)) + + +def ask(candidates: dict[str, list[str]], layer_order: list[str], + default_net: str, e1_label: str, e2_label: str, + contact1: str, contact2: str, + buildup_layers: list[str] | None = None) -> Selection | None: + """Show the dialog; returns None on cancel.""" + app = QApplication.instance() or QApplication([]) + dlg = _Dialog(candidates, layer_order, default_net, e1_label, e2_label, + contact1, contact2, buildup_layers or []) + dlg.raise_() + dlg.activateWindow() + if dlg.exec() != QDialog.Accepted: + return None + return dlg.selection() diff --git a/fill_resistance/errors.py b/fill_resistance/errors.py new file mode 100644 index 0000000..3b4dac8 --- /dev/null +++ b/fill_resistance/errors.py @@ -0,0 +1,34 @@ +"""User-facing error hierarchy. + +Every UserFacingError message is shown both on stdout (KiCad status bar) +and in a matplotlib error figure, so keep messages self-contained and +actionable. +""" + + +class UserFacingError(Exception): + pass + + +class ApiVersionError(UserFacingError): + pass + + +class SelectionError(UserFacingError): + pass + + +class CandidateError(UserFacingError): + pass + + +class ElectrodeError(UserFacingError): + pass + + +class ConnectivityError(UserFacingError): + pass + + +class GridSizeError(UserFacingError): + pass diff --git a/fill_resistance/geometry.py b/fill_resistance/geometry.py new file mode 100644 index 0000000..8e167ac --- /dev/null +++ b/fill_resistance/geometry.py @@ -0,0 +1,317 @@ +"""Plain geometry data model. No kipy imports here. + +Everything is int64 nanometers in KiCad board coordinates (y grows down); +z grows from the board top surface downwards through the stackup. +Problem is the complete solver input and doubles as the JSON dump schema, +so the whole pipeline downstream of board_io runs without KiCad. + +Schema v2 is multi-layer: per-layer fills at stackup depths, linked by +via/through-pad barrels. v1 dumps (single layer, no vias) still load. +""" +from __future__ import annotations + +import json +import math +from dataclasses import dataclass, field +from pathlib import Path + +import numpy as np + +JSON_SCHEMA_VERSION = 4 + + +@dataclass(frozen=True) +class Rect: + x0: int + y0: int + x1: int + y1: int + layer_name: str + + @classmethod + def normalized(cls, xa: int, ya: int, xb: int, yb: int, layer_name: str) -> "Rect": + return cls(min(xa, xb), min(ya, yb), max(xa, xb), max(ya, yb), layer_name) + + @property + def width(self) -> int: + return self.x1 - self.x0 + + @property + def height(self) -> int: + return self.y1 - self.y0 + + +@dataclass +class Polygon: + outline: np.ndarray # (N, 2) int64 nm, open ring + holes: list[np.ndarray] = field(default_factory=list) + + +@dataclass +class LayerFill: + layer_name: str + thickness_nm: int + z_nm: int # copper center depth from board top + polygons: list[Polygon] + + +@dataclass +class SurfaceBuildup: + """Solder (plus optional added copper) sitting on an outer copper + layer inside solder-mask openings (zones on F.Mask/B.Mask).""" + layer_name: str # copper layer it sits on + polygons: list[Polygon] + + +@dataclass +class Electrode: + """One PART of a current-injection terminal: a drawn rectangle or a + selected pad. A terminal (V+ or V-) is a LIST of parts, all merged + into one equipotential contact (externally bonded). `polygons` + (board nm) is the exact copper shape when known (pads); None means + the rectangle itself is the shape. `contact` = 'all' or a layer + name: which included layers this part touches.""" + rect: Rect # bounding box (labels/summary) + contact: str = "all" + polygons: list[Polygon] | None = None + label: str = "rect" + + +@dataclass +class ViaLink: + """A conductive barrel (via or plated through-hole pad) linking copper + layers whose z lies within [z_top_nm, z_bot_nm].""" + x: int + y: int + drill_nm: int + z_top_nm: int + z_bot_nm: int + kind: str = "via" # "via" | "pad" + + def spans(self, z_nm: int) -> bool: + return self.z_top_nm - 1 <= z_nm <= self.z_bot_nm + 1 + + def barrel_resistance(self, length_nm: int, rho_ohm_m: float, + plating_nm: int) -> float: + """Barrel segment resistance over length_nm: thin-wall annulus of + plating around the drill.""" + area_m2 = math.pi * (self.drill_nm * 1e-9) * (plating_nm * 1e-9) + return rho_ohm_m * (length_nm * 1e-9) / area_m2 + + +@dataclass +class Problem: + board_path: str + net_name: str + rho_ohm_m: float + plating_nm: int + layers: list[LayerFill] # sorted by z_nm (top first) + vias: list[ViaLink] + electrodes1: list[Electrode] # V+ terminal parts (merged) + electrodes2: list[Electrode] # V- terminal parts (merged) + thickness_source: str = "stackup" + buildups: list[SurfaceBuildup] = field(default_factory=list) + solder_thickness_nm: int = 50_000 + solder_rho_ohm_m: float = 1.32e-7 + extra_cu_nm: int = 0 + + @property + def layer_names(self) -> list[str]: + return [l.layer_name for l in self.layers] + + def sigma_s(self, layer_index: int) -> float: + """Sheet conductance of one layer [S per square].""" + return (self.layers[layer_index].thickness_nm * 1e-9) / self.rho_ohm_m + + def copper_bbox(self) -> tuple[int, int, int, int]: + xs = np.concatenate([p.outline[:, 0] + for l in self.layers for p in l.polygons]) + ys = np.concatenate([p.outline[:, 1] + for l in self.layers for p in l.polygons]) + return int(xs.min()), int(ys.min()), int(xs.max()), int(ys.max()) + + +def arc_points(start, mid, end, tol_nm: float) -> np.ndarray: + """Tessellate a start/mid/end arc into points from start (inclusive) + to end (exclusive), max sagitta <= tol_nm. Collinear input degrades + to just the start point (straight segment).""" + sx, sy = float(start[0]), float(start[1]) + mx, my = float(mid[0]), float(mid[1]) + ex, ey = float(end[0]), float(end[1]) + + d = 2.0 * (sx * (my - ey) + mx * (ey - sy) + ex * (sy - my)) + chord = math.hypot(ex - sx, ey - sy) + if abs(d) < 1e-9 * max(chord, 1.0): + return np.array([[start[0], start[1]]], dtype=np.int64) + ux = ((sx**2 + sy**2) * (my - ey) + (mx**2 + my**2) * (ey - sy) + + (ex**2 + ey**2) * (sy - my)) / d + uy = ((sx**2 + sy**2) * (ex - mx) + (mx**2 + my**2) * (sx - ex) + + (ex**2 + ey**2) * (mx - sx)) / d + r = math.hypot(sx - ux, sy - uy) + + a0 = math.atan2(sy - uy, sx - ux) + a1 = math.atan2(my - uy, mx - ux) + a2 = math.atan2(ey - uy, ex - ux) + two_pi = 2.0 * math.pi + d01 = (a1 - a0) % two_pi + d02 = (a2 - a0) % two_pi + sweep = d02 if d01 <= d02 else d02 - two_pi + + tol = min(tol_nm, 0.999 * r) + dtheta_max = 2.0 * math.acos(1.0 - tol / r) + n = max(2, int(math.ceil(abs(sweep) / dtheta_max))) + ks = np.arange(n) + angs = a0 + sweep * ks / n + pts = np.stack([ux + r * np.cos(angs), uy + r * np.sin(angs)], axis=1) + return np.round(pts).astype(np.int64) + + +def linearize_ring(nodes: list, tol_nm: float) -> np.ndarray: + """nodes: list of ('pt', (x, y)) or ('arc', (start, mid, end)) tuples, + already in board nm. Returns an (N, 2) int64 open ring.""" + parts = [] + for kind, data in nodes: + if kind == "pt": + parts.append(np.array([[data[0], data[1]]], dtype=np.int64)) + elif kind == "arc": + parts.append(arc_points(data[0], data[1], data[2], tol_nm)) + else: + raise ValueError(f"unknown polyline node kind: {kind}") + ring = np.concatenate(parts, axis=0) + if len(ring) > 1 and (ring[0] == ring[-1]).all(): + ring = ring[:-1] + return ring + + +# --- JSON dump / load ------------------------------------------------------- + +def _poly_to_json(p: Polygon) -> dict: + return {"outline": p.outline.tolist(), "holes": [h.tolist() for h in p.holes]} + + +def _poly_from_json(d: dict) -> Polygon: + return Polygon(outline=np.asarray(d["outline"], dtype=np.int64), + holes=[np.asarray(h, dtype=np.int64) for h in d["holes"]]) + + +def _electrode_to_json(e: Electrode) -> dict: + return { + "rect": vars(e.rect) | {}, + "contact": e.contact, + "label": e.label, + "polygons": (None if e.polygons is None + else [_poly_to_json(poly) for poly in e.polygons]), + } + + +def _electrode_from_json(d: dict) -> Electrode: + return Electrode( + rect=_rect_from_json(d["rect"]), + contact=d.get("contact", "all"), + label=d.get("label", "rect"), + polygons=(None if d.get("polygons") is None + else [_poly_from_json(pd) for pd in d["polygons"]]), + ) + + +def problem_to_json(p: Problem) -> dict: + return { + "schema_version": JSON_SCHEMA_VERSION, + "board_path": p.board_path, + "net_name": p.net_name, + "rho_ohm_m": p.rho_ohm_m, + "plating_nm": p.plating_nm, + "thickness_source": p.thickness_source, + "electrodes1": [_electrode_to_json(e) for e in p.electrodes1], + "electrodes2": [_electrode_to_json(e) for e in p.electrodes2], + "layers": [ + { + "layer_name": l.layer_name, + "thickness_nm": l.thickness_nm, + "z_nm": l.z_nm, + "polygons": [_poly_to_json(poly) for poly in l.polygons], + } + for l in p.layers + ], + "vias": [vars(v) | {} for v in p.vias], + "buildups": [ + {"layer_name": b.layer_name, + "polygons": [_poly_to_json(poly) for poly in b.polygons]} + for b in p.buildups + ], + "solder_thickness_nm": p.solder_thickness_nm, + "solder_rho_ohm_m": p.solder_rho_ohm_m, + "extra_cu_nm": p.extra_cu_nm, + } + + +def _rect_from_json(rd: dict) -> Rect: + return Rect(int(rd["x0"]), int(rd["y0"]), int(rd["x1"]), int(rd["y1"]), + rd["layer_name"]) + + +def problem_from_json(d: dict) -> Problem: + version = d.get("schema_version", 1) + if version == 1: + # v1: single layer, no vias, rect electrodes + return Problem( + board_path=d["board_path"], + net_name=d["net_name"], + rho_ohm_m=float(d["rho_ohm_m"]), + plating_nm=18_000, + layers=[LayerFill( + layer_name=d["layer_name"], + thickness_nm=int(d["thickness_nm"]), + z_nm=0, + polygons=[_poly_from_json(pd) for pd in d["polygons"]], + )], + vias=[], + electrodes1=[Electrode(rect=_rect_from_json(d["rect1"]))], + electrodes2=[Electrode(rect=_rect_from_json(d["rect2"]))], + thickness_source=d.get("thickness_source", "unknown"), + ) + return Problem( + board_path=d["board_path"], + net_name=d["net_name"], + rho_ohm_m=float(d["rho_ohm_m"]), + plating_nm=int(d["plating_nm"]), + layers=[ + LayerFill( + layer_name=ld["layer_name"], + thickness_nm=int(ld["thickness_nm"]), + z_nm=int(ld["z_nm"]), + polygons=[_poly_from_json(pd) for pd in ld["polygons"]], + ) + for ld in d["layers"] + ], + vias=[ + ViaLink(x=int(vd["x"]), y=int(vd["y"]), drill_nm=int(vd["drill_nm"]), + z_top_nm=int(vd["z_top_nm"]), z_bot_nm=int(vd["z_bot_nm"]), + kind=vd.get("kind", "via")) + for vd in d["vias"] + ], + electrodes1=( + [_electrode_from_json(ed) for ed in d["electrodes1"]] + if version >= 3 else [_electrode_from_json(d["electrode1"])]), + electrodes2=( + [_electrode_from_json(ed) for ed in d["electrodes2"]] + if version >= 3 else [_electrode_from_json(d["electrode2"])]), + thickness_source=d.get("thickness_source", "unknown"), + buildups=[ + SurfaceBuildup( + layer_name=bd["layer_name"], + polygons=[_poly_from_json(pd) for pd in bd["polygons"]]) + for bd in d.get("buildups", []) + ], + solder_thickness_nm=int(d.get("solder_thickness_nm", 50_000)), + solder_rho_ohm_m=float(d.get("solder_rho_ohm_m", 1.32e-7)), + extra_cu_nm=int(d.get("extra_cu_nm", 0)), + ) + + +def save_problem(p: Problem, path: Path) -> None: + path.write_text(json.dumps(problem_to_json(p)), encoding="utf-8") + + +def load_problem(path: Path) -> Problem: + return problem_from_json(json.loads(Path(path).read_text(encoding="utf-8"))) diff --git a/fill_resistance/main.py b/fill_resistance/main.py new file mode 100644 index 0000000..1f99c18 --- /dev/null +++ b/fill_resistance/main.py @@ -0,0 +1,108 @@ +"""Top-level orchestration for the KiCad-launched action. + +Flow: connect -> read the two selected contacts (rectangles/pads) -> +gather fills -> selection dialog (net, layers, contacts, current, cell) +-> extract vias -> solve -> figures + report. + +Every failure is reported twice: on stdout (lands in the KiCad status-bar +warning list) and as a matplotlib error figure, so it cannot be missed. +""" +from __future__ import annotations + +import sys +import traceback + +from . import config, pipeline, report +from .errors import CandidateError, UserFacingError + + +def _fail(message: str, outdir) -> None: + print(f"ERROR: {message}") + from . import plots + fig = plots.fig_error(message) + plots.save_and_show([(fig, "error")], outdir) + sys.exit(1) + + +def main() -> None: + outdir = None + try: + from kipy.errors import ApiError + + from . import board_io, dialog + try: + kicad, board = board_io.connect() + stackup = board_io.get_stackup_info(board) + es1, es2, net_hint = board_io.get_electrodes(board) + if board_io.any_zone_unfilled(board) or config.ALWAYS_REFILL: + board_io.refill(board) + fills = board_io.gather_net_fills(board) + candidate_nets = board_io.nets_overlapping(fills, es1, es2) + buildups = board_io.gather_mask_buildups(board) + except ApiError as e: + raise UserFacingError( + f"KiCad API error: {e}\nIf KiCad is showing a dialog, close " + f"it and run again." + ) + + if not candidate_nets: + raise CandidateError( + "No copper zone fill overlaps both contacts. Check that both " + "sit over (or in) filled pours and that the fills are up to " + "date (press B in the board editor)." + ) + + def group_label(parts): + names = [p.label for p in parts[:3]] + more = f" +{len(parts) - 3}" if len(parts) > 3 else "" + return f"{len(parts)}× " + ", ".join(names) + more + + def group_contact(parts): + contacts = {p.contact for p in parts} + return contacts.pop() if len(contacts) == 1 else "auto" + + default_net = (net_hint if net_hint in candidate_nets + else candidate_nets[0]) + selection = dialog.ask( + candidates={n: list(fills[n].keys()) for n in candidate_nets}, + layer_order=stackup.names, + default_net=default_net, + e1_label=group_label(es1), e2_label=group_label(es2), + contact1=group_contact(es1), contact2=group_contact(es2), + buildup_layers=sorted(buildups.keys()), + ) + if selection is None: + print("cancelled") + return + + if selection.contact1 != "auto": + for e in es1: + e.contact = selection.contact1 + if selection.contact2 != "auto": + for e in es2: + e.contact = selection.contact2 + if selection.cell_um is not None: + config.CELL_UM_OVERRIDE = selection.cell_um + + try: + problem = board_io.build_problem( + board, selection.net, selection.layers, es1, es2, stackup, + fills, + buildups=(buildups if selection.include_buildup else None), + extra_cu_um=selection.extra_cu_um) + outdir = report.make_output_dir(board_io.board_dir(board)) + except ApiError as e: + raise UserFacingError(f"KiCad API error: {e}") + + report.write_geometry_dump(outdir, problem) + pipeline.run(problem, outdir, show=True, i_test=selection.current_a, + freq_hz=selection.freq_hz, + contact_model=selection.contact_model) + except UserFacingError as e: + _fail(str(e), outdir) + except Exception: + _fail(traceback.format_exc(), outdir) + + +if __name__ == "__main__": + main() diff --git a/fill_resistance/pipeline.py b/fill_resistance/pipeline.py new file mode 100644 index 0000000..ac52f25 --- /dev/null +++ b/fill_resistance/pipeline.py @@ -0,0 +1,51 @@ +"""Geometry-in -> results-out pipeline shared by the KiCad entrypoint and +the offline standalone runner.""" +from __future__ import annotations + +from pathlib import Path + +from . import config, plots, raster, report, solver +from .geometry import Problem +from .solver import Result + + +def run(problem: Problem, outdir: Path | None, show: bool = True, + i_test: float | None = None, freq_hz: float = 0.0, + contact_model: str | None = None) -> Result: + if i_test is None: + i_test = config.TEST_CURRENT_A + h = raster.choose_cell_size(problem.copper_bbox(), len(problem.layers)) + print(f"rasterizing {len(problem.layers)} layer(s) at cell size " + f"{h / 1000:.1f} um ...") + stack = raster.rasterize_stack(problem, h) + print(f"grid {stack.shape2d[1]}x{stack.shape2d[0]}x{stack.nlayers}, " + f"{int(stack.masks.sum())} copper cells, {len(problem.vias)} " + f"via/pad barrel(s)") + + e1, e2 = raster.electrode_masks(stack, problem) + parts1, parts2 = raster.electrode_partition(stack, problem) + + print(f"solving @ {i_test:g} A" + + (f", {freq_hz:g} Hz" if freq_hz > 0 else " DC") + " ...") + result = solver.run_solve(problem, stack, e1, e2, i_test, freq_hz, + contact_model, parts1, parts2) + for prefix, pcs in (("P", result.part_currents1), + ("N", result.part_currents2)): + for i, (label, amps) in enumerate(pcs): + print(f" {prefix}{i + 1} ({label}): {amps:.4g} A " + f"({100 * amps / i_test:.1f}%)") + + if outdir is not None: + outdir.mkdir(parents=True, exist_ok=True) + report.write_summary(outdir, problem, stack, result) + print(report.result_line(result, problem, stack)) + + figs = [ + (plots.fig_raster(stack, e1, e2, problem, result), "1_raster_map"), + (plots.fig_potential(result, stack, e1, e2, problem), "2_potential"), + (plots.fig_current(result, stack, e1, e2, problem), + "3_current_density"), + (plots.fig_power(result, stack, e1, e2, problem), "4_power_density"), + ] + plots.save_and_show(figs, outdir, show=show) + return result diff --git a/fill_resistance/plots.py b/fill_resistance/plots.py new file mode 100644 index 0000000..e66dc63 --- /dev/null +++ b/fill_resistance/plots.py @@ -0,0 +1,391 @@ +"""Figures: per-layer rasterized maps, potential, current density, power +density, and the error figure. PNGs are saved BEFORE any window opens. + +Backend: interactive if a GUI toolkit exists (tkinter, else Qt), else Agg +with os.startfile on the saved PNGs so results are never silent. +""" +from __future__ import annotations + +import os +import textwrap +from pathlib import Path + +import matplotlib +import numpy as np + + +def _pick_backend(): + """matplotlib.use() is lazy and 'succeeds' for backends whose GUI + toolkit is missing (KiCad's Python has no tkinter), so probe the + toolkits explicitly.""" + try: + import tkinter # noqa: F401 + return "TkAgg" + except Exception: + pass + for qt in ("PySide6", "PyQt6", "PyQt5", "PySide2"): + try: + __import__(qt) + return "QtAgg" if qt in ("PySide6", "PyQt6") else "Qt5Agg" + except Exception: + continue + return None + + +INTERACTIVE_BACKEND = _pick_backend() +matplotlib.use(INTERACTIVE_BACKEND or "Agg") + +import matplotlib.pyplot as plt # noqa: E402 (after backend selection) +from matplotlib.colors import ListedColormap, LogNorm # noqa: E402 +from matplotlib.patches import Patch # noqa: E402 + +from . import config # noqa: E402 + +_BG = "#f5f3f0" +_COPPER = "#c98b4e" +_E1_COLOR = "#c8385a" +_E2_COLOR = "#2f6fb0" +_VIA_COLOR = "#2d6b45" +_SOLDER = "#9aa3ad" # tin-gray: solder buildup areas +_INK = "#3a3a3a" +_GRID_INK = "#b8b4ae" + + +def _fmt_si(value: float, unit: str) -> str: + for scale, prefix in ((1.0, ""), (1e-3, "m"), (1e-6, "µ")): + if abs(value) >= scale: + return f"{value / scale:.4g} {prefix}{unit}" + return f"{value:.3g} {unit}" + + +def _suptitle(problem, stack, result=None) -> str: + ny, nx = stack.shape2d + parts = [] + if result is not None: + parts.append(f"R = {result.R_ohm * 1000:.4g} mΩ") + parts.append(f"P = {_fmt_si(result.P_total, 'W')} @ " + f"{result.i_test:g} A") + if result.freq_hz > 0: + parts.append(f"f = {result.freq_hz / 1e3:g} kHz " + f"(δ={result.skin_depth_um:.0f} µm, lower bound)") + parts.append(problem.net_name) + parts.append(f"{nx}×{ny}×{stack.nlayers} @ {stack.h_nm / 1000:.0f} µm") + return " | ".join(parts) + + +def _layer_fig(stack, window_title: str): + L = stack.nlayers + ny, nx = stack.shape2d + aspect = ny / nx + w = 9.5 + row_h = min(max(w * aspect * 0.9 + 0.6, 1.8), 8.5 / L) + fig, axes = plt.subplots(L, 1, figsize=(w, row_h * L + 1.4), + sharex=True, sharey=True, squeeze=False) + axes = axes[:, 0] + if INTERACTIVE_BACKEND: + fig.canvas.manager.set_window_title(window_title) + for ax, name in zip(axes, stack.layer_names): + ax.set_ylabel(f"{name}\ny [mm]", fontsize=8) + ax.tick_params(colors=_INK, labelsize=8) + for s in ax.spines.values(): + s.set_color(_GRID_INK) + axes[-1].set_xlabel("x [mm]") + return fig, axes + + +def _electrode_labels(ax, stack, e1_l, e2_l): + """Label each connected contact part (multi-part terminals get one + label per island, largest first, up to 4).""" + from scipy import ndimage + for e, label, color in ((e1_l, "V+", _E1_COLOR), (e2_l, "V−", _E2_COLOR)): + if not e.any(): + continue + labels, n = ndimage.label(e) + sizes = ndimage.sum_labels(np.ones_like(labels), labels, + range(1, n + 1)) + order = np.argsort(sizes)[::-1][:4] + 1 + for comp in order: + ii, jj = np.nonzero(labels == comp) + cx = (stack.x0_nm + (jj.mean() + 0.5) * stack.h_nm) * 1e-6 + cy = (stack.y0_nm + (ii.mean() + 0.5) * stack.h_nm) * 1e-6 + ax.annotate(label, (cx, cy), xytext=(0, 0), + textcoords="offset points", color="white", + fontsize=9, fontweight="bold", ha="center", + va="center", + bbox=dict(boxstyle="round,pad=0.2", fc=color, + ec="none", alpha=0.9)) + + +def _via_markers(ax, problem, layer): + xs = [v.x * 1e-6 for v in problem.vias if v.spans(layer.z_nm)] + ys = [v.y * 1e-6 for v in problem.vias if v.spans(layer.z_nm)] + if xs: + ax.plot(xs, ys, ".", ms=2.5, color=_VIA_COLOR, alpha=0.7) + + +def area_tag(sign: str, index: int) -> str: + """Short injection-area tag: P1, P2, ... for V+; N1, N2, ... for V-.""" + return f"{'P' if sign == '+' else 'N'}{index + 1}" + + +def _injection_area_labels(ax, li, layer_name, problem, result): + """Mark every injection area with its short tag (currents live in + the legend).""" + groups = ((problem.electrodes1, "+", _E1_COLOR), + (problem.electrodes2, "-", _E2_COLOR)) + for parts, sign, color in groups: + for i, el in enumerate(parts): + if el.contact != "all" and el.contact != layer_name: + continue + cx = (el.rect.x0 + el.rect.x1) / 2e6 + cy = (el.rect.y0 + el.rect.y1) / 2e6 + ax.annotate(area_tag(sign, i), (cx, cy), xytext=(0, 0), + textcoords="offset points", color="white", + fontsize=8, fontweight="bold", ha="center", + va="center", + bbox=dict(boxstyle="round,pad=0.15", fc=color, + ec="none", alpha=0.9)) + + +def fig_raster(stack, e1, e2, problem, result=None): + fig, axes = _layer_fig(stack, "Fill Resistance - rasterized map") + cmap = ListedColormap([_BG, _COPPER, _E1_COLOR, _E2_COLOR, _SOLDER]) + has_buildup = stack.buildup is not None and stack.buildup.any() + for li, ax in enumerate(axes): + codes = np.zeros(stack.shape2d, dtype=np.uint8) + codes[stack.masks[li]] = 1 + if has_buildup: + codes[stack.buildup[li]] = 4 + codes[e1[li]] = 2 + codes[e2[li]] = 3 + ax.imshow(codes, cmap=cmap, vmin=0, vmax=4, origin="upper", + extent=stack.extent_mm(), interpolation="nearest") + _via_markers(ax, problem, problem.layers[li]) + if result is not None and (result.part_currents1 + or result.part_currents2): + _injection_area_labels(ax, li, stack.layer_names[li], problem, + result) + else: + _electrode_labels(ax, stack, e1[li], e2[li]) + handles = [Patch(fc=_COPPER, label="copper"), + Patch(fc=_VIA_COLOR, label="vias")] + if has_buildup: + handles.append(Patch( + fc=_SOLDER, + label=f"solder buildup " + f"({problem.solder_thickness_nm / 1000:.0f} µm" + + (f" + {problem.extra_cu_nm / 1000:.0f} µm Cu" + if problem.extra_cu_nm else "") + ")")) + if result is not None and (result.part_currents1 + or result.part_currents2): + entries = ([("+", _E1_COLOR, i, amps) + for i, (_, amps) in enumerate(result.part_currents1)] + + [("-", _E2_COLOR, i, amps) + for i, (_, amps) in enumerate(result.part_currents2)]) + shown = entries[:14] + for sign, color, i, amps in shown: + handles.append(Patch( + fc=color, + label=f"{area_tag(sign, i)}: {amps:.3g} A " + f"({100 * amps / result.i_test:.0f}%)")) + if len(entries) > len(shown): + handles.append(Patch(fc="#00000000", + label=f"... +{len(entries) - len(shown)} " + f"more in summary.txt")) + else: + handles += [Patch(fc=_E1_COLOR, label="V+"), + Patch(fc=_E2_COLOR, label="V−")] + axes[0].legend(handles=handles, loc="upper right", fontsize=7, + framealpha=0.9) + fig.suptitle("Rasterized fill + electrodes | " + + _suptitle(problem, stack, result), fontsize=10, color=_INK) + fig.tight_layout() + return fig + + +def fig_potential(result, stack, e1, e2, problem): + fig, axes = _layer_fig(stack, "Fill Resistance - potential") + vmax = float(np.nanmax(result.V)) + unit, scale = ("mV", 1e3) if vmax < 0.1 else ("V", 1.0) + cmap = matplotlib.colormaps[config.CMAP_POTENTIAL].copy() + cmap.set_bad(_BG) + im = None + for li, ax in enumerate(axes): + vs = result.V[li] * scale + im = ax.imshow(vs, cmap=cmap, vmin=0, vmax=vmax * scale, + origin="upper", extent=stack.extent_mm(), + interpolation="nearest") + if np.isfinite(vs).sum() > 4: + ext = stack.extent_mm() + ny, nx = stack.shape2d + xs = np.linspace(ext[0], ext[1], nx, endpoint=False) + xs += (xs[1] - xs[0]) / 2 + ys = np.linspace(ext[3], ext[2], ny, endpoint=False) + ys += (ys[1] - ys[0]) / 2 + with np.errstate(invalid="ignore"): + ax.contour(xs, ys, vs, levels=15, colors="white", + linewidths=0.4, alpha=0.5) + _electrode_labels(ax, stack, e1[li], e2[li]) + cb = fig.colorbar(im, ax=axes, shrink=0.85) + cb.set_label(f"potential [{unit}] @ {result.i_test:g} A", fontsize=9) + fig.suptitle("Potential | " + _suptitle(problem, stack, result), + fontsize=10, color=_INK) + return fig + + +def _field_fig(result, stack, e1, e2, problem, data3, cmap_name, dyn_range, + label, title, window): + """Shared per-layer LogNorm field figure (current, power).""" + fig, axes = _layer_fig(stack, window) + vmax = float(np.nanmax(data3)) + cmap = matplotlib.colormaps[cmap_name].copy() + cmap.set_bad(_BG) + if config.LOG_CURRENT_SCALE and vmax > 0: + norm = LogNorm(vmin=vmax / dyn_range, vmax=vmax) + else: + norm = None + im = None + for li, ax in enumerate(axes): + d = data3[li] + shown = np.clip(d, vmax / dyn_range, None) if norm is not None else d + im = ax.imshow(shown, cmap=cmap, norm=norm, origin="upper", + extent=stack.extent_mm(), interpolation="nearest") + _electrode_labels(ax, stack, e1[li], e2[li]) + if vmax > 0: + li, i, j = np.unravel_index(np.nanargmax(data3), data3.shape) + mx = (stack.x0_nm + (j + 0.5) * stack.h_nm) * 1e-6 + my = (stack.y0_nm + (i + 0.5) * stack.h_nm) * 1e-6 + axes[li].plot(mx, my, "o", ms=9, mfc="none", mec="white", mew=1.4) + axes[li].annotate(f"max {vmax:.3g}", (mx, my), xytext=(10, -10), + textcoords="offset points", color="white", + fontsize=8, + bbox=dict(boxstyle="round,pad=0.2", fc="#00000088", + ec="none")) + cb = fig.colorbar(im, ax=axes, shrink=0.85) + cb.set_label(label, fontsize=9) + fig.suptitle(title + " | " + _suptitle(problem, stack, result), + fontsize=10, color=_INK) + return fig, axes + + +def fig_current(result, stack, e1, e2, problem): + fig, axes = _field_fig( + result, stack, e1, e2, problem, result.Jmag * 1e-6, + config.CMAP_CURRENT, config.CURRENT_DYNAMIC_RANGE, + f"|J| [A/mm²] @ {result.i_test:g} A", + "Current density (log)", "Fill Resistance - current density") + # mark the hottest via + if result.via_reports: + v = result.via_reports[0] + for ax in axes: + ax.plot(v.x_mm, v.y_mm, "s", ms=7, mfc="none", mec="#7fe0a8", + mew=1.2) + axes[0].annotate( + f"hottest via {v.current_a:.3g} A", (v.x_mm, v.y_mm), + xytext=(10, 10), textcoords="offset points", color="white", + fontsize=8, + bbox=dict(boxstyle="round,pad=0.2", fc="#2d6b45", ec="none")) + return fig + + +def fig_power(result, stack, e1, e2, problem): + # W/m^2 -> W/mm^2 + fig, axes = _field_fig( + result, stack, e1, e2, problem, result.Parea * 1e-6, + config.CMAP_POWER, config.POWER_DYNAMIC_RANGE, + f"p [W/mm²] @ {result.i_test:g} A", + "Power density (log)", "Fill Resistance - power density") + for li, ax in enumerate(axes): + ax.set_title(f"P({stack.layer_names[li]}) = " + f"{_fmt_si(result.P_layers[li], 'W')}", + fontsize=8, color=_INK, loc="right", pad=2) + return fig + + +def fig_error(message: str): + fig, ax = plt.subplots(figsize=(9, 4.5)) + ax.axis("off") + ax.set_title("Fill Resistance — ERROR", color="#b02a2a", + fontsize=14, fontweight="bold", loc="left") + wrapped = "\n".join( + textwrap.fill(line, width=90) for line in message.splitlines() + ) + ax.text(0.0, 0.95, wrapped, family="monospace", fontsize=9, + va="top", ha="left", color=_INK, transform=ax.transAxes) + fig.tight_layout() + return fig + + +def _resolve_label_overlaps(fig): + """Measure every annotation's rendered box and greedily push + overlapping labels upward until nothing collides. Runs on the real + renderer, so it handles any font/DPI.""" + from matplotlib.text import Annotation + try: + fig.canvas.draw() + renderer = fig.canvas.get_renderer() + except Exception: + return + for ax in fig.axes: + anns = [c for c in ax.get_children() if isinstance(c, Annotation)] + placed = [] + for a in sorted(anns, key=lambda t: t.get_window_extent(renderer).x0): + try: + bb = a.get_window_extent(renderer) + except Exception: + continue + guard = 50 + while guard > 0: + hit = next((p for p in placed if bb.overlaps(p)), None) + if hit is None: + break + push_px = (hit.y1 - bb.y0) + 3.0 + dx, dy = a.xyann + a.xyann = (dx, dy + push_px * 72.0 / fig.dpi) + bb = a.get_window_extent(renderer) + guard -= 1 + placed.append(bb) + + +def _raise_windows(): + """Best effort: bring plot windows in front of KiCad (windows spawned + by a background process tend to open behind).""" + for num in plt.get_fignums(): + try: + win = plt.figure(num).canvas.manager.window + if hasattr(win, "attributes"): # Tk + win.attributes("-topmost", True) + win.after(300, lambda w=win: w.attributes("-topmost", False)) + else: # Qt + win.raise_() + win.activateWindow() + except Exception: + pass + + +def save_and_show(figs_named: list[tuple], outdir: Path | None, + show: bool = True) -> list[Path]: + """figs_named: [(figure, basename), ...]. Saves first, then shows.""" + saved = [] + for fig, _ in figs_named: + _resolve_label_overlaps(fig) + if outdir is not None: + outdir.mkdir(parents=True, exist_ok=True) + for fig, name in figs_named: + p = outdir / f"{name}.png" + fig.savefig(p, dpi=config.DPI, facecolor="white", + bbox_inches="tight") + saved.append(p) + print(f"saved {p}") + if show and config.INTERACTIVE: + if INTERACTIVE_BACKEND: + _raise_windows() + plt.show() + else: + for p in saved: + try: + os.startfile(p) # windows: open in default viewer + except Exception: + pass + plt.close("all") + return saved diff --git a/fill_resistance/raster.py b/fill_resistance/raster.py new file mode 100644 index 0000000..a0904aa --- /dev/null +++ b/fill_resistance/raster.py @@ -0,0 +1,263 @@ +"""Rasterization of the fill polygons onto a shared multi-layer grid, and +electrode mask construction. + +Grid convention: layer l, row i, col j maps to the cell center + x = x0_nm + (j + 0.5) * h_nm + y = y0_nm + (i + 0.5) * h_nm +in KiCad board coordinates (y grows down). Row 0 is the minimum-y row, +the TOP of the board as drawn in the editor; plots use origin='upper'. +All layers share the same frame, so cell (i, j) is vertically aligned +across layers (via links connect equal (i, j) on different layers). + +Connectivity restriction lives in solver.py: it needs the via edges. +""" +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np +from matplotlib.path import Path as MplPath +from scipy import ndimage + +from . import config +from .errors import ElectrodeError, GridSizeError +from .geometry import Electrode, Problem, Rect + +# 4-connectivity: matches the in-plane 5-point stencil of the solver +_STRUCT4 = ndimage.generate_binary_structure(2, 1) + + +@dataclass +class RasterStack: + masks: np.ndarray # bool (L, ny, nx), True = copper + x0_nm: float # grid origin (outer corner of cell [., 0, 0]) + y0_nm: float + h_nm: float + layer_names: list[str] + buildup: np.ndarray | None = None # bool (L, ny, nx): solder buildup + # (mask opening ∩ copper) + + @property + def nlayers(self) -> int: + return self.masks.shape[0] + + @property + def shape2d(self) -> tuple[int, int]: + return self.masks.shape[1:] + + def cell_centers(self, i0: int, i1: int, j0: int, j1: int): + xs = self.x0_nm + (np.arange(j0, j1) + 0.5) * self.h_nm + ys = self.y0_nm + (np.arange(i0, i1) + 0.5) * self.h_nm + return np.meshgrid(xs, ys) + + def cell_of(self, x_nm: float, y_nm: float) -> tuple[int, int] | None: + """(i, j) of the cell containing the point, or None if outside.""" + ny, nx = self.shape2d + j = int((x_nm - self.x0_nm) / self.h_nm) + i = int((y_nm - self.y0_nm) / self.h_nm) + if 0 <= i < ny and 0 <= j < nx: + return i, j + return None + + def extent_mm(self) -> tuple[float, float, float, float]: + """imshow extent (left, right, bottom, top) for origin='upper', + y axis in board orientation (increasing downward).""" + ny, nx = self.shape2d + return ( + self.x0_nm * 1e-6, + (self.x0_nm + nx * self.h_nm) * 1e-6, + (self.y0_nm + ny * self.h_nm) * 1e-6, + self.y0_nm * 1e-6, + ) + + +def choose_cell_size(bbox_nm: tuple[int, int, int, int], nlayers: int) -> float: + """Pick the cell size h [nm]; TARGET_CELLS counts TOTAL cells across + all layers. Raise if the grid would exceed HARD_MAX_CELLS.""" + x0, y0, x1, y1 = bbox_nm + w, ht = float(x1 - x0), float(y1 - y0) + if w <= 0 or ht <= 0: + raise GridSizeError("Copper geometry has a degenerate bounding box.") + + if config.CELL_UM_OVERRIDE is not None: + h = config.CELL_UM_OVERRIDE * 1000.0 + else: + h = math.sqrt(w * ht * nlayers / config.TARGET_CELLS) + h = min(max(h, config.MIN_CELL_UM * 1000.0), config.MAX_CELL_UM * 1000.0) + + ncells = math.ceil(w / h) * math.ceil(ht / h) * nlayers + if ncells > config.HARD_MAX_CELLS: + raise GridSizeError( + f"Grid would need ~{ncells / 1e6:.1f} M cells over {nlayers} " + f"layer(s) at cell size {h / 1000:.0f} um (limit " + f"{config.HARD_MAX_CELLS / 1e6:.0f} M). Raise MAX_CELL_UM / " + f"CELL_UM_OVERRIDE in config.py, deselect layers, or measure a " + f"smaller region." + ) + return h + + +def _paint_ring(stack: RasterStack, ring: np.ndarray, value: bool, + target: np.ndarray) -> None: + """Set target (2D) cells whose center lies inside ring to `value`, + testing only cells within the ring's bbox (cheap for small holes).""" + ny, nx = stack.shape2d + h = stack.h_nm + j0 = max(0, int((ring[:, 0].min() - stack.x0_nm) / h) - 1) + j1 = min(nx, int((ring[:, 0].max() - stack.x0_nm) / h) + 2) + i0 = max(0, int((ring[:, 1].min() - stack.y0_nm) / h) - 1) + i1 = min(ny, int((ring[:, 1].max() - stack.y0_nm) / h) + 2) + if i0 >= i1 or j0 >= j1: + return + xg, yg = stack.cell_centers(i0, i1, j0, j1) + pts = np.column_stack([xg.ravel(), yg.ravel()]) + # Path(closed=True) treats the LAST vertex as the CLOSEPOLY dummy, so + # the first vertex must be appended or the ring loses its last corner + verts = np.vstack([ring, ring[:1]]) + inside = MplPath(verts, closed=True).contains_points(pts) + inside = inside.reshape(i1 - i0, j1 - j0) + sub = target[i0:i1, j0:j1] + sub[inside] = value + + +def rasterize_stack(problem: Problem, h_nm: float) -> RasterStack: + """Rasterize every included layer onto one shared frame.""" + x0, y0, x1, y1 = problem.copper_bbox() + m = config.MARGIN_CELLS + nx = math.ceil((x1 - x0) / h_nm) + 2 * m + ny = math.ceil((y1 - y0) / h_nm) + 2 * m + stack = RasterStack( + masks=np.zeros((len(problem.layers), ny, nx), dtype=bool), + x0_nm=x0 - m * h_nm, + y0_nm=y0 - m * h_nm, + h_nm=h_nm, + layer_names=problem.layer_names, + ) + for li, layer in enumerate(problem.layers): + for poly in layer.polygons: + pmask = np.zeros((ny, nx), dtype=bool) + _paint_ring(stack, poly.outline, True, pmask) + for hole in poly.holes: + _paint_ring(stack, hole, False, pmask) + stack.masks[li] |= pmask + + if problem.buildups: + stack.buildup = np.zeros_like(stack.masks) + index = {name: li for li, name in enumerate(stack.layer_names)} + for b in problem.buildups: + li = index.get(b.layer_name) + if li is None: + continue + for poly in b.polygons: + pmask = np.zeros((ny, nx), dtype=bool) + _paint_ring(stack, poly.outline, True, pmask) + for hole in poly.holes: + _paint_ring(stack, hole, False, pmask) + stack.buildup[li] |= pmask + stack.buildup &= stack.masks # solder wets exposed copper only + return stack + + +def _rect_cells(stack: RasterStack, rect: Rect) -> np.ndarray: + """Bool (ny, nx) mask of cells whose center lies inside the rectangle.""" + ny, nx = stack.shape2d + h = stack.h_nm + out = np.zeros((ny, nx), dtype=bool) + j0 = max(0, int(math.ceil((rect.x0 - stack.x0_nm) / h - 0.5))) + j1 = min(nx, int(math.floor((rect.x1 - stack.x0_nm) / h - 0.5)) + 1) + i0 = max(0, int(math.ceil((rect.y0 - stack.y0_nm) / h - 0.5))) + i1 = min(ny, int(math.floor((rect.y1 - stack.y0_nm) / h - 0.5)) + 1) + if i0 < i1 and j0 < j1: + out[i0:i1, j0:j1] = True + return out + + +def _electrode_cells2d(stack: RasterStack, e: Electrode) -> np.ndarray: + """2D footprint of the electrode shape (pad polygons or rectangle).""" + if e.polygons: + cells = np.zeros(stack.shape2d, dtype=bool) + for poly in e.polygons: + pm = np.zeros(stack.shape2d, dtype=bool) + _paint_ring(stack, poly.outline, True, pm) + for hole in poly.holes: + _paint_ring(stack, hole, False, pm) + cells |= pm + if not cells.any(): + # shape smaller than one grid cell (small pad): use the cell + # containing its center + r = e.rect + c = stack.cell_of((r.x0 + r.x1) / 2, (r.y0 + r.y1) / 2) + if c is not None: + cells[c] = True + return cells + return _rect_cells(stack, e.rect) + + +def electrode_masks(stack: RasterStack, problem: Problem + ) -> tuple[np.ndarray, np.ndarray]: + """Terminal mask = OR over its parts; part = shape ∩ copper on the + part's contact layer(s). contact 'all' = every included layer (bolted + lug / through pad); a layer name = that layer only. Every part must + individually land on copper (clear feedback). V+/V- must not overlap; + touching is checked later, only for the equipotential contact model.""" + def build(parts: list[Electrode], which: str) -> np.ndarray: + e = np.zeros_like(stack.masks) + for el in parts: + cells2d = _electrode_cells2d(stack, el) + part = np.zeros_like(stack.masks) + for li, name in enumerate(stack.layer_names): + if el.contact == "all" or el.contact == name: + part[li] = cells2d & stack.masks[li] + if not part.any(): + raise ElectrodeError( + f"A {which} contact part ({el.label}) does not overlap " + f"any copper of the selected fill on contact layer(s) " + f"'{el.contact}' (or is smaller than one grid cell)." + ) + e |= part + if not e.any(): + raise ElectrodeError(f"The {which} terminal has no contact parts.") + return e + + e1 = build(problem.electrodes1, "V+") + e2 = build(problem.electrodes2, "V-") + + if (e1 & e2).any(): + raise ElectrodeError( + "The V+ and V- contact areas overlap on the copper grid. " + "Move them apart." + ) + return e1, e2 + + +def electrode_partition(stack: RasterStack, problem: Problem + ) -> tuple[list, list]: + """Per-part cell masks for both terminals, as [(label, mask3d), ...]. + Cells covered by several overlapping parts are attributed to the + FIRST part (first-wins partition), so part currents sum exactly to + the terminal current.""" + def build(parts: list[Electrode]) -> list: + out = [] + claimed = np.zeros_like(stack.masks) + for el in parts: + cells2d = _electrode_cells2d(stack, el) + m = np.zeros_like(stack.masks) + for li, name in enumerate(stack.layer_names): + if el.contact == "all" or el.contact == name: + m[li] = cells2d & stack.masks[li] + m &= ~claimed + claimed |= m + out.append((el.label, m)) + return out + + return build(problem.electrodes1), build(problem.electrodes2) + + +def electrodes_touch(stack: RasterStack, e1: np.ndarray, + e2: np.ndarray) -> str | None: + """Layer name where the terminals are 4-adjacent, or None.""" + for li in range(stack.nlayers): + if (ndimage.binary_dilation(e1[li], structure=_STRUCT4) & e2[li]).any(): + return stack.layer_names[li] + return None diff --git a/fill_resistance/report.py b/fill_resistance/report.py new file mode 100644 index 0000000..f186e16 --- /dev/null +++ b/fill_resistance/report.py @@ -0,0 +1,146 @@ +"""Output directory, summary.txt, geometry dump, stdout one-liner.""" +from __future__ import annotations + +from datetime import datetime +from pathlib import Path + +import numpy as np + +from . import config +from .geometry import Problem, save_problem +from .raster import RasterStack +from .solver import Result + + +def make_output_dir(board_dir: Path) -> Path: + stamp = datetime.now().strftime("%Y%m%d-%H%M%S") + out = Path(board_dir) / config.OUTPUT_DIRNAME / stamp + out.mkdir(parents=True, exist_ok=True) + return out + + +def write_geometry_dump(outdir: Path, problem: Problem) -> Path: + p = outdir / "geometry_dump.json" + save_problem(problem, p) + return p + + +def result_line(result: Result, problem: Problem, stack: RasterStack) -> str: + ny, nx = stack.shape2d + ac = (f" @ {result.freq_hz / 1e3:g} kHz (lower bound)" + if result.freq_hz > 0 else "") + return (f"R = {result.R_ohm * 1000:.4g} mOhm{ac}, " + f"P = {result.P_total:.4g} W @ {result.i_test:g} A " + f"(net {problem.net_name}, {'+'.join(stack.layer_names)}, " + f"grid {nx}x{ny}x{stack.nlayers}, cell {stack.h_nm / 1000:.0f} um)") + + +def _electrode_line(e) -> str: + r = e.rect + return (f"{e.label:12s} contact={e.contact:8s} " + f"x [{r.x0 / 1e6:.2f}, {r.x1 / 1e6:.2f}] " + f"y [{r.y0 / 1e6:.2f}, {r.y1 / 1e6:.2f}] mm") + + +def write_summary(outdir: Path, problem: Problem, stack: RasterStack, + result: Result) -> Path: + ny, nx = stack.shape2d + info = result.solve_info + lines = [ + "fill_resistance summary", + "=======================", + f"board: {problem.board_path}", + f"net: {problem.net_name}", + f"test current: {result.i_test:g} A", + f"resistivity: {problem.rho_ohm_m:.3e} ohm*m", + f"via plating: {problem.plating_nm / 1000:.0f} um", + "", + (f"frequency: " + + (f"{result.freq_hz:g} Hz (skin depth {result.skin_depth_um:.0f} um)" + if result.freq_hz > 0 else "DC")), + f"RESISTANCE: {result.R_ohm * 1000:.6g} mOhm" + + (" (AC LOWER BOUND: lateral/proximity redistribution not modeled)" + if result.freq_hz > 0 else ""), + f"VOLTAGE DROP: {result.R_ohm * result.i_test * 1000:.4g} mV " + f"@ {result.i_test:g} A", + f"TOTAL POWER: {result.P_total:.6g} W @ {result.i_test:g} A", + f" in vias: {result.P_vias:.4g} W", + f" power balance: {result.power_balance_rel:.2e} (consistency)", + "", + "layers (top to bottom):", + ] + if problem.buildups and stack.buildup is not None: + eq_um = (problem.solder_thickness_nm / 1000 + * problem.rho_ohm_m / problem.solder_rho_ohm_m + + problem.extra_cu_nm / 1000) + cell_mm2 = (stack.h_nm * 1e-6) ** 2 + per_layer = {name: float(stack.buildup[li].sum()) * cell_mm2 + for li, name in enumerate(stack.layer_names) + if stack.buildup[li].any()} + areas = ", ".join(f"{n}: {a:.0f} mm^2" for n, a in per_layer.items()) + lines.insert(-1, f"solder buildup: " + f"{problem.solder_thickness_nm / 1000:.0f} um solder" + + (f" + {problem.extra_cu_nm / 1000:.0f} um Cu" + if problem.extra_cu_nm else "") + + f" = {eq_um:.1f} um equivalent Cu ({areas})") + for li, layer in enumerate(problem.layers): + ac = (f" Rs_AC/Rs_DC={result.rs_ratios[li]:.2f}" + if result.freq_hz > 0 else "") + lines.append( + f" {layer.layer_name:8s} t={layer.thickness_nm / 1000:5.1f} um " + f"z={layer.z_nm / 1000:7.1f} um " + f"P={result.P_layers[li]:.4g} W " + f"maxJ={float(np.nanmax(result.Jmag[li])) * 1e-6 if np.isfinite(result.Jmag[li]).any() else 0:.4g} A/mm^2" + + ac + ) + lines += [ + "", + f"grid: {nx} x {ny} x {stack.nlayers} cells @ " + f"{stack.h_nm / 1000:.1f} um", + f"copper cells: {int(stack.masks.sum())}", + f"free unknowns: {result.n_free}", + f"solver: {info.method}" + + (f", {info.iterations} iters, residual {info.residual:.2e}" + if info.iterations is not None else ""), + f"I1/I2 @ 1V: {result.I1_a:.9g} / {result.I2_a:.9g} A " + f"(mismatch {result.mismatch_rel:.2e})", + f"timings [s]: " + f"{', '.join(f'{k}={v:.2f}' for k, v in result.timings.items())}", + "", + f"contact model: {result.contact_model}" + + (" (uniform orthogonal injection; R is the upper contact bound)" + if result.contact_model == "uniform" else " (ideal bonded lug)"), + f"terminals:", + f" V+ ({len(problem.electrodes1)} injection area(s)):", + *(f" {_electrode_line(e)}" for e in problem.electrodes1), + f" V- ({len(problem.electrodes2)} injection area(s)):", + *(f" {_electrode_line(e)}" for e in problem.electrodes2), + ] + if result.part_currents1 or result.part_currents2: + how = ("prescribed by area share (uniform model)" + if result.contact_model == "uniform" + else "computed flux (equipotential model)") + lines += ["", f"current per injection area @ {result.i_test:g} A " + f"({how}):"] + for sign, pcs in (("+", result.part_currents1), + ("-", result.part_currents2)): + for i, (label, amps) in enumerate(pcs): + tag = f"{'P' if sign == '+' else 'N'}{i + 1}" + lines.append(f" {tag:4s} {label:24s} {amps:9.4g} A " + f"({100 * amps / result.i_test:5.1f}%)") + if result.via_reports: + n_shown = min(10, len(result.via_reports)) + lines += [ + "", + f"vias/pads carrying current (top {n_shown} of " + f"{len(result.via_reports)}, @ {result.i_test:g} A):", + " x [mm] y [mm] kind drill I [A] P [W]", + ] + for v in result.via_reports[:n_shown]: + lines.append( + f" {v.x_mm:8.2f} {v.y_mm:8.2f} {v.kind:5s} " + f"{v.drill_mm:5.2f} {v.current_a:8.4g} {v.power_w:.4g}" + ) + p = outdir / "summary.txt" + p.write_text("\n".join(lines), encoding="utf-8") + return p diff --git a/fill_resistance/skin.py b/fill_resistance/skin.py new file mode 100644 index 0000000..e27a332 --- /dev/null +++ b/fill_resistance/skin.py @@ -0,0 +1,78 @@ +"""Skin-effect corrections: frequency-dependent effective sheet +resistance of a copper foil and via-barrel wall. + +1D diffusion through the foil thickness (exact): with tau = (1+j)/delta, +the internal impedance per square of a foil of thickness t is + + one-sided field (plane over a return plane): Zs = tau*rho * coth(tau*t) + two-sided field (isolated foil): Zs = tau*rho/2 * coth(tau*t/2) + +Both reduce to rho/t at DC and to rho/delta (resp. rho/(2*delta)) at +high frequency. R_AC = Re(Zs) is used as the effective sheet resistance. + +HONESTY NOTE (also in the README): only the through-thickness current +crowding is modeled. Lateral redistribution (proximity effect - AC +current following the minimum-inductance path) needs a magneto- +quasistatic solve and is NOT captured; since the resistance-driven +distribution is the minimum-dissipation one, the reported AC resistance +is a rigorous LOWER BOUND at the given frequency. +""" +from __future__ import annotations + +import cmath +import math + +MU0 = 4e-7 * math.pi + + +def skin_depth_m(freq_hz: float, rho_ohm_m: float) -> float: + return math.sqrt(2.0 * rho_ohm_m / (2.0 * math.pi * freq_hz * MU0)) + + +def _coth(x: complex) -> complex: + return 1.0 / cmath.tanh(x) + + +def sheet_resistance_ac(thickness_m: float, freq_hz: float, + rho_ohm_m: float, sides: int = 1) -> float: + """Effective sheet resistance [ohm/sq] of a foil at freq_hz. + sides=1: field on one side (plane facing a return plane, conservative); + sides=2: symmetric field on both sides (isolated foil).""" + if freq_hz <= 0.0: + return rho_ohm_m / thickness_m + delta = skin_depth_m(freq_hz, rho_ohm_m) + tau = (1.0 + 1.0j) / delta + if sides == 2: + zs = tau * rho_ohm_m / 2.0 * _coth(tau * thickness_m / 2.0) + else: + zs = tau * rho_ohm_m * _coth(tau * thickness_m) + return zs.real + + +def resistance_factor(thickness_m: float, freq_hz: float, + rho_ohm_m: float, sides: int = 1) -> float: + """R_AC / R_DC of a foil (or barrel wall) of the given thickness.""" + if freq_hz <= 0.0: + return 1.0 + return (sheet_resistance_ac(thickness_m, freq_hz, rho_ohm_m, sides) + / (rho_ohm_m / thickness_m)) + + +def parse_frequency(text: str) -> float: + """'0', '100k', '1.5M', '142500' -> Hz. Empty/invalid -> 0 (DC).""" + t = text.strip().lower().replace(",", ".").removesuffix("hz").strip() + if not t: + return 0.0 + mult = 1.0 + if t.endswith("meg"): + mult, t = 1e6, t[:-3] + elif t.endswith("m"): + mult, t = 1e6, t[:-1] + elif t.endswith("k"): + mult, t = 1e3, t[:-1] + elif t.endswith("g"): + mult, t = 1e9, t[:-1] + try: + return max(0.0, float(t) * mult) + except ValueError: + return 0.0 diff --git a/fill_resistance/solver.py b/fill_resistance/solver.py new file mode 100644 index 0000000..8129da7 --- /dev/null +++ b/fill_resistance/solver.py @@ -0,0 +1,549 @@ +"""Coupled multi-layer finite-difference solver. + +Each included copper layer is a 2D 5-point sheet with per-layer face +conductance sigma_s = t/rho [S] (square cells: independent of h); via and +plated-through-pad barrels add vertical conductances between vertically +aligned cells of the layers they span AND reach copper on. A barrel +passing an antipad still bridges the layers above/below it with the full +barrel length. At freq > 0 the per-layer sheet conductances and the +barrel walls get the 1D skin-effect correction (see skin.py; AC results +are a rigorous lower bound - lateral redistribution is not modeled). + +Two contact models for the terminals (each terminal = merged parts): + +- "uniform" (default): a conductor pressed onto the contact area injects + the current orthogonally with UNIFORM surface density: every contact + cell sources (sinks) I/N. The in-plane current density ramps across + the contact instead of being zero. The pure-Neumann system is grounded + at one V- cell (that cell's sink share is exactly the flux that exits + through the ground reference, so the solution equals the singular + system's). R = ( - ) / I; because + the injection and averaging weights coincide, sum(edge powers) = I^2 R + holds exactly and remains the consistency check. + +- "equipotential": ideal bonded lug; contact cells are Dirichlet + (V+ = 1 V, V- = 0). R from the exact discrete electrode flux. Touching + terminals are rejected (a direct face would short the Dirichlet + regions); with "uniform" contacts touching is physically fine. + +The two models bracket a real contact: R_equipotential <= R_real <= +R_uniform. Missing neighbors give no matrix term = insulated boundary. +Current density per layer comes from face currents (np.gradient across +the NaN staircase boundary would pollute the field). Power density per +layer distributes each in-plane edge's dissipation half to each endpoint +cell. All reported fields are rescaled to the test current I_test. +""" +from __future__ import annotations + +import time +from dataclasses import dataclass, field + +import numpy as np +from scipy import sparse +from scipy.sparse import csgraph +from scipy.sparse import linalg as sla + +from . import config, skin +from .errors import ConnectivityError, ElectrodeError +from .geometry import Problem +from .raster import RasterStack, electrodes_touch + + +@dataclass +class SolveInfo: + method: str # "spsolve" | "cg+jacobi" + n_unknowns: int + iterations: int | None = None + residual: float | None = None + + +@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 + + +@dataclass +class ViaReport: + x_mm: float + y_mm: float + kind: str + drill_mm: float + current_a: float # max barrel-segment current @ I_test + power_w: float # total barrel dissipation @ I_test + + +@dataclass +class Result: + R_ohm: float + i_test: float + V: np.ndarray # (L, ny, nx) volts @ I_test, NaN off-copper + Jmag: np.ndarray # (L, ny, nx) A/m^2 @ I_test + Parea: np.ndarray # (L, ny, nx) W/m^2 @ I_test + layer_names: list[str] + P_total: float # I_test^2 * R + P_layers: list[float] # in-plane dissipation per layer @ I_test + P_vias: float # total barrel dissipation @ I_test + power_balance_rel: float # |sum(edge powers) - I^2 R| / I^2 R + via_reports: list[ViaReport] # sorted by current, descending + I1_a: float # electrode currents (unit drive) + I2_a: float + mismatch_rel: float + n_free: int + solve_info: SolveInfo + # per-part terminal currents @ I_test: [(label, amps), ...]; + # computed flux for "equipotential", prescribed area share for "uniform" + part_currents1: list = field(default_factory=list) + part_currents2: list = field(default_factory=list) + contact_model: str = "uniform" + freq_hz: float = 0.0 + 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) + + +def _shifts2d(): + return [ + ((slice(None), slice(None, -1)), (slice(None), slice(1, None))), + ((slice(None, -1), slice(None)), (slice(1, None), slice(None))), + ] + + +def _sigma_2d(stack: RasterStack, li: int, sigma_layer: float, + sigma_buildup: float) -> np.ndarray | None: + """Per-cell sheet conductance for one layer, or None if uniform.""" + if stack.buildup is None or sigma_buildup <= 0 \ + or not stack.buildup[li].any(): + return None + s = np.full(stack.shape2d, sigma_layer) + s[stack.buildup[li]] += sigma_buildup + return s + + +def build_edges(stack: RasterStack, problem: Problem, sigmas: list[float], + via_factor: float = 1.0, + sigma_buildup: float = 0.0) -> Edges: + """All copper-copper conductances: in-plane faces + via barrels. + sigmas: effective (possibly AC) sheet conductance per layer; + via_factor: R_AC/R_DC of the barrel wall; sigma_buildup: extra sheet + conductance on solder-buildup cells. Faces between cells of unequal + conductance use the harmonic mean (series half-cells), which reduces + exactly to sigma for uniform regions.""" + L, ny, nx = stack.masks.shape + plane = ny * nx + aa, bb, ww, vv = [], [], [], [] + + for li in range(L): + m = stack.masks[li] + sig = sigmas[li] + scell = _sigma_2d(stack, li, sig, sigma_buildup) + base = li * plane + for src, dst in _shifts2d(): + pair = m[src] & m[dst] + ii, jj = np.nonzero(pair) + if src[0] == slice(None): # horizontal: j, j+1 + a = base + ii * nx + jj + b = a + 1 + else: # vertical: i, i+1 + a = base + ii * nx + jj + b = a + nx + aa.append(a.astype(np.int64)) + bb.append(b.astype(np.int64)) + if scell is None: + ww.append(np.full(len(a), sig)) + else: + s_a = scell[src][pair] + s_b = scell[dst][pair] + ww.append(2.0 * s_a * s_b / (s_a + s_b)) + vv.append(np.full(len(a), -1, dtype=np.int32)) + + for vi, via in enumerate(problem.vias): + cell = stack.cell_of(via.x, via.y) + if cell is None: + continue + i, j = cell + present = [li for li, layer in enumerate(problem.layers) + if via.spans(layer.z_nm) and stack.masks[li, i, j]] + for la, lb 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 + i * nx + j], dtype=np.int64)) + bb.append(np.array([lb * plane + i * nx + j], dtype=np.int64)) + ww.append(np.array([1.0 / r])) + vv.append(np.array([vi], dtype=np.int32)) + + if not aa: + raise ConnectivityError("No copper found on the selected layers.") + return Edges(a=np.concatenate(aa), b=np.concatenate(bb), + w=np.concatenate(ww), via_index=np.concatenate(vv)) + + +def connected_restrict(stack: RasterStack, e1: np.ndarray, e2: np.ndarray, + edges: Edges) -> bool: + """Keep only components (through-plane AND through-via) touching both + terminals. Mutates stack.masks / e1 / e2. Returns True if anything + was dropped (caller must rebuild edges).""" + 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) + common = np.intersect1d(np.unique(labels3[e1]), np.unique(labels3[e2])) + 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." + ) + keep = np.isin(labels3, common) & stack.masks + changed = bool((stack.masks & ~keep).any()) + stack.masks &= keep + e1 &= keep + e2 &= keep + return changed + + +def _assemble(state: np.ndarray, edges: Edges, rhs_extra: np.ndarray | 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.""" + 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." + ) + + free = state == 1 + n_free = int(free.sum()) + if n_free == 0: + raise ElectrodeError( + "No free copper cells remain between the terminals - the " + "contacts cover the whole fill at this grid resolution." + ) + idx = np.full(n, -1, dtype=np.int64) + idx[free] = np.arange(n_free) + + diag = np.zeros(n_free) + rhs = np.zeros(n_free) + 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 rhs_extra is not None: + rhs += rhs_extra[free] + + ff = fa & fb + rows = np.concatenate([idx[edges.a[ff]], idx[edges.b[ff]], + np.arange(n_free)]) + cols = np.concatenate([idx[edges.b[ff]], idx[edges.a[ff]], + np.arange(n_free)]) + vals = np.concatenate([-edges.w[ff], -edges.w[ff], diag]) + A = sparse.coo_matrix((vals, (rows, cols)), + shape=(n_free, n_free)).tocsr() + return A, rhs, idx + + +def solve_system(A: sparse.csr_matrix, b: np.ndarray) -> tuple[np.ndarray, SolveInfo]: + n = A.shape[0] + if n <= config.SPSOLVE_MAX_UNKNOWNS: + x = sla.spsolve(A.tocsc(), b) + return x, SolveInfo(method="spsolve", n_unknowns=n) + + # The matrix is SPD, so CG is guaranteed to converge. Jacobi is the + # only preconditioner in scipy that keeps the preconditioned operator + # SPD without a factorization that can break down at this scale. + d = A.diagonal() + M = sla.LinearOperator((n, n), lambda v: v / d) + iters = 0 + + def count(_): + nonlocal iters + iters += 1 + + try: + x, code = sla.cg(A, b, M=M, rtol=config.CG_TOL, + maxiter=config.CG_MAXITER, callback=count) + except TypeError: # scipy < 1.12 uses tol= + x, code = sla.cg(A, b, M=M, tol=config.CG_TOL, + maxiter=config.CG_MAXITER, callback=count) + if code != 0: + raise RuntimeError( + f"CG did not converge in {config.CG_MAXITER} iterations " + f"(code {code}). Try a coarser grid or raise CG_MAXITER." + ) + res = float(np.linalg.norm(b - A @ x) / np.linalg.norm(b)) + return x, SolveInfo(method="cg+jacobi", n_unknowns=n, iterations=iters, + residual=res) + + +def _face_current_density(V2: np.ndarray, mask2: np.ndarray, sigma: float, + h_m: float, t_m: float, + sig2d: np.ndarray | None = None, + rho: float | None = None) -> np.ndarray: + """|J| (A/m^2) for one layer from face currents; V2 in volts. + With a per-cell conductance map (buildup), face currents use the + harmonic mean and J is referenced to the conductance-equivalent + copper thickness t_eq = sigma_cell * rho (equals the geometric t for + plain DC copper).""" + ny, nx = mask2.shape + face_x = mask2[:, :-1] & mask2[:, 1:] + face_y = mask2[:-1, :] & mask2[1:, :] + if sig2d is None: + wx = wy = sigma + teq = np.full((ny, nx), t_m) + else: + wx = 2.0 * sig2d[:, :-1] * sig2d[:, 1:] / (sig2d[:, :-1] + sig2d[:, 1:]) + wy = 2.0 * sig2d[:-1, :] * sig2d[1:, :] / (sig2d[:-1, :] + sig2d[1:, :]) + teq = sig2d * rho + with np.errstate(invalid="ignore"): + Ix = np.where(face_x, (V2[:, :-1] - V2[:, 1:]) * wx, 0.0) + Iy = np.where(face_y, (V2[:-1, :] - V2[1:, :]) * wy, 0.0) + IxP = np.zeros((ny, nx + 1)) + IxP[:, 1:nx] = Ix + IyP = np.zeros((ny + 1, nx)) + IyP[1:ny, :] = Iy + Jx = 0.5 * (IxP[:, :-1] + IxP[:, 1:]) + Jy = 0.5 * (IyP[:-1, :] + IyP[1:, :]) + Jmag = np.hypot(Jx, Jy) / (h_m * teq) + Jmag[~mask2] = np.nan + 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) + x, info = solve_system(A, rhs) + + Vflat = np.zeros(state.size) + Vflat[state == 2] = 1.0 + Vflat[state == 1] = x + + Ie = edges.w * (Vflat[edges.a] - Vflat[edges.b]) + sa, sb = state[edges.a], state[edges.b] + I1 = float(Ie[sa == 2].sum() - Ie[sb == 2].sum()) + I2 = float(Ie[sb == 3].sum() - Ie[sa == 3].sum()) + mismatch = abs(I1 - I2) / max(abs(I1), abs(I2), 1e-300) + R = 1.0 / (0.5 * (I1 + I2)) + 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) + x, info = solve_system(A, rhs) + + Vflat = np.zeros(n) + Vflat[state == 1] = x + + v_plus = float(Vflat[e1f].mean()) + v_minus = float(Vflat[e2f].mean()) + R = (v_plus - v_minus) / 1.0 + Vflat = Vflat - v_minus # display reference: = 0 + + # quality: KCL residual of the solved system + res = info.residual + if res is None: + res = float(np.linalg.norm(A @ x - rhs) + / max(np.linalg.norm(rhs), 1e-300)) + return Vflat, R, 1.0, 1.0, res, 1.0, info + + +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 + discrete flux out of the part's cells (same-terminal internal edges + carry zero, opposite-terminal edges are forbidden). Uniform: the + injection is prescribed, so a part carries exactly its cell share.""" + out = [] + for label, mask3 in parts: + pf = mask3.ravel() & e_flat + n = int(pf.sum()) + if contact_model == "uniform": + amps = i_test * n / max(n_terminal_cells, 1) + else: + ina = pf[edges.a] + inb = pf[edges.b] + amps = abs(float(Ie[ina].sum() - Ie[inb].sum())) * scale + out.append((label, amps)) + 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 + sigmas = [ + 1.0 / skin.sheet_resistance_ac( + problem.layers[li].thickness_nm * 1e-9, freq_hz, + problem.rho_ohm_m, config.SKIN_SIDES) + for li in range(L) + ] + rs_ratios = [ + skin.resistance_factor(problem.layers[li].thickness_nm * 1e-9, + freq_hz, problem.rho_ohm_m, config.SKIN_SIDES) + for li in range(L) + ] + via_factor = skin.resistance_factor(problem.plating_nm * 1e-9, freq_hz, + problem.rho_ohm_m, sides=2) + sigma_buildup = 0.0 + if problem.buildups and stack.buildup is not None \ + and stack.buildup.any(): + sigma_buildup = 1.0 / skin.sheet_resistance_ac( + problem.solder_thickness_nm * 1e-9, freq_hz, + problem.solder_rho_ohm_m, config.SKIN_SIDES) + if problem.extra_cu_nm > 0: + sigma_buildup += 1.0 / skin.sheet_resistance_ac( + problem.extra_cu_nm * 1e-9, freq_hz, problem.rho_ohm_m, + config.SKIN_SIDES) + eq_um = sigma_buildup * problem.rho_ohm_m * 1e6 + print(f"solder buildup: {problem.solder_thickness_nm / 1000:.0f} um " + f"solder + {problem.extra_cu_nm / 1000:.0f} um Cu on " + f"{int(stack.buildup.sum())} cells " + f"(= {eq_um:.1f} um equivalent copper)") + if freq_hz > 0: + depth = skin.skin_depth_m(freq_hz, problem.rho_ohm_m) + print(f"AC @ {freq_hz:g} Hz: skin depth {depth * 1e6:.0f} um, " + f"per-layer Rs ratio " + f"{', '.join(f'{r:.2f}' for r in rs_ratios)}, " + f"via factor {via_factor:.2f}") + + t0 = time.perf_counter() + edges = build_edges(stack, problem, sigmas, via_factor, sigma_buildup) + if connected_restrict(stack, e1, e2, edges): + edges = build_edges(stack, problem, sigmas, via_factor, sigma_buildup) + if stack.buildup is not None: + stack.buildup &= stack.masks + for _, m in (parts1 or []) + (parts2 or []): + m &= stack.masks # follow the component restriction + timings["edges_s"] = time.perf_counter() - t0 + + t0 = time.perf_counter() + if contact_model == "equipotential": + Vflat, R, I1, I2, mismatch, volts_per_amp, info = \ + _solve_equipotential(stack, e1, e2, edges) + else: + Vflat, R, I1, I2, mismatch, volts_per_amp, info = \ + _solve_uniform(stack, e1, e2, edges) + timings["solve_s"] = time.perf_counter() - t0 + + 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()) + P_total = i_test ** 2 * R + balance = abs((sum(P_layers) + P_vias) - P_total) / max(P_total, 1e-300) + + # 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, + contact_model, int(e1.sum())) + 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 + V3 = np.full((L, ny, nx), np.nan) + V3[stack.masks] = Vflat.reshape(L, ny, nx)[stack.masks] * s + J3 = np.stack([ + _face_current_density( + np.nan_to_num(V3[li]), stack.masks[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) + ]) + timings["postprocess_s"] = time.perf_counter() - t0 + + return 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/standalone.py b/fill_resistance/standalone.py new file mode 100644 index 0000000..3785cca --- /dev/null +++ b/fill_resistance/standalone.py @@ -0,0 +1,80 @@ +"""Offline runner: solve a geometry_dump.json without KiCad. + + python -m fill_resistance.standalone dump.json [--current 40] + [--cell-um 50] [--layers F.Cu,In1.Cu] [--no-show] [--out DIR] + [--force-iterative] + +This is the dev loop and the convergence-study tool (KiCad 10 has no +headless API server, so the plugin path always needs the GUI). +""" +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +from . import config, pipeline +from .errors import UserFacingError +from .geometry import load_problem + + +def main(argv=None) -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("dump", type=Path, help="geometry_dump.json from a plugin run") + ap.add_argument("--current", type=float, default=None, + help="test current [A] (default: config TEST_CURRENT_A)") + ap.add_argument("--freq", type=str, default="0", + help="frequency, e.g. 142k or 1.5M (default: DC). " + "AC results are a lower bound (skin per foil only)") + ap.add_argument("--cell-um", type=float, default=None, + help="force grid cell size [um]") + ap.add_argument("--layers", type=str, default=None, + help="comma-separated subset of layers to include") + ap.add_argument("--out", type=Path, default=None, + help="output directory (default: next to the dump)") + ap.add_argument("--no-show", action="store_true", + help="save PNGs only, no windows") + ap.add_argument("--contact-model", choices=["uniform", "equipotential"], + default=None, help="contact model (default: config)") + ap.add_argument("--strip-buildup", action="store_true", + help="ignore solder buildup stored in the dump") + ap.add_argument("--extra-cu-um", type=float, default=None, + help="override the added copper in mask openings [um]") + ap.add_argument("--force-iterative", action="store_true", + help="use CG (Jacobi) regardless of problem size") + args = ap.parse_args(argv) + + if args.cell_um is not None: + config.CELL_UM_OVERRIDE = args.cell_um + if args.no_show: + config.INTERACTIVE = False + if args.force_iterative: + config.SPSOLVE_MAX_UNKNOWNS = 0 + + problem = load_problem(args.dump) + if args.strip_buildup: + problem.buildups = [] + if args.extra_cu_um is not None: + problem.extra_cu_nm = int(args.extra_cu_um * 1000) + if args.layers: + keep = [s.strip() for s in args.layers.split(",")] + problem.layers = [l for l in problem.layers if l.layer_name in keep] + if not problem.layers: + print(f"ERROR: no layer of the dump matches --layers {args.layers}", + file=sys.stderr) + return 1 + + from .skin import parse_frequency + outdir = args.out if args.out is not None else args.dump.parent + try: + pipeline.run(problem, outdir, show=not args.no_show, + i_test=args.current, freq_hz=parse_frequency(args.freq), + contact_model=args.contact_model) + except UserFacingError as e: + print(f"ERROR: {e}", file=sys.stderr) + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/icons/fill_res_24_dark.png b/icons/fill_res_24_dark.png new file mode 100644 index 0000000..e0ec6b1 Binary files /dev/null and b/icons/fill_res_24_dark.png differ diff --git a/icons/fill_res_24_light.png b/icons/fill_res_24_light.png new file mode 100644 index 0000000..4f7c1d0 Binary files /dev/null and b/icons/fill_res_24_light.png differ diff --git a/plugin.json b/plugin.json new file mode 100644 index 0000000..93a2e99 --- /dev/null +++ b/plugin.json @@ -0,0 +1,21 @@ +{ + "$schema": "https://go.kicad.org/api/schemas/v1", + "identifier": "th.co.b4l.fill-resistance", + "name": "Fill Resistance", + "description": "DC resistance of a copper zone fill between two rectangle electrodes", + "runtime": { + "type": "python" + }, + "actions": [ + { + "identifier": "fill-resistance-run", + "name": "Fill Resistance", + "description": "Select two rectangles marking the contact areas, then run to compute the fill resistance between them", + "entrypoint": "fill_res_action.py", + "show-button": true, + "scopes": ["pcb"], + "icons-light": ["icons/fill_res_24_light.png"], + "icons-dark": ["icons/fill_res_24_dark.png"] + } + ] +} diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..919a38e --- /dev/null +++ b/requirements.txt @@ -0,0 +1,5 @@ +kicad-python>=0.7.0 +numpy +scipy +matplotlib +PySide6 diff --git a/smoke/smoke_probe.py b/smoke/smoke_probe.py new file mode 100644 index 0000000..a195197 --- /dev/null +++ b/smoke/smoke_probe.py @@ -0,0 +1,100 @@ +"""M1 smoke probe: verify the IPC API surface against a live KiCad. + +Run from the dev venv while KiCad is open with the board loaded: + .venv\\Scripts\\python.exe smoke\\smoke_probe.py + +Deliberately uses kipy directly (not board_io) so it works even if +board_io has a bug. Prints: version, board path, selection contents, +shapes summary, zones + filled polygon stats, stackup table. +""" +import sys +import traceback + +from kipy import KiCad +from kipy.board_types import BoardRectangle, BoardShape +from kipy.proto.board.board_pb2 import BoardStackupLayerType +from kipy.proto.board.board_types_pb2 import ZoneType +from kipy.util.board_layer import canonical_name + + +def mm(nm): + return nm / 1e6 + + +def main(): + print("== connect ==") + kicad = KiCad() + kicad.ping() + print("version:", kicad.get_version()) + try: + print("check_version:", kicad.check_version()) + except Exception as e: + print("check_version raised:", e) + + board = kicad.get_board() + print("board.name:", board.name) + print("board_filename:", getattr(board.document, "board_filename", "?")) + + print("\n== selection ==") + sel = list(board.get_selection()) + print(f"{len(sel)} item(s) selected") + for item in sel: + line = f" {type(item).__name__}" + if isinstance(item, BoardRectangle): + tl, br = item.top_left, item.bottom_right + line += (f" layer={canonical_name(item.layer)}" + f" tl=({mm(tl.x):.2f}, {mm(tl.y):.2f})mm" + f" br=({mm(br.x):.2f}, {mm(br.y):.2f})mm") + elif isinstance(item, BoardShape): + line += f" layer={canonical_name(item.layer)}" + print(line) + + print("\n== shapes (board-wide rectangles) ==") + shapes = list(board.get_shapes()) + rect_shapes = [s for s in shapes if isinstance(s, BoardRectangle)] + print(f"{len(shapes)} shapes total, {len(rect_shapes)} rectangles") + for s in rect_shapes[:20]: + tl, br = s.top_left, s.bottom_right + print(f" rect on {canonical_name(s.layer)}: " + f"({mm(tl.x):.2f}, {mm(tl.y):.2f}) - ({mm(br.x):.2f}, {mm(br.y):.2f}) mm") + + print("\n== zones ==") + for zone in board.get_zones(): + ztype = ZoneType.Name(zone.type) + net = zone.net.name if zone.net is not None else "" + layers = [canonical_name(l) for l in zone.layers] + print(f" zone '{zone.name}' type={ztype} net={net} " + f"layers={layers} filled={zone.filled}") + try: + for layer, polys in zone.filled_polygons.items(): + narcs = 0 + nnodes = 0 + nholes = 0 + for p in polys: + nnodes += len(p.outline.nodes) + nholes += len(p.holes) + narcs += sum(1 for n in p.outline.nodes if n.has_arc) + print(f" fill on {canonical_name(layer)}: {len(polys)} " + f"poly(s), {nnodes} outline nodes, {nholes} holes, " + f"{narcs} arc nodes") + except Exception: + print(" filled_polygons FAILED:") + traceback.print_exc() + + print("\n== stackup ==") + try: + for sl in board.get_stackup().layers: + tname = BoardStackupLayerType.Name(sl.type) + lname = canonical_name(sl.layer) if sl.type == \ + BoardStackupLayerType.BSLT_COPPER else "-" + print(f" {tname:18s} layer={lname:8s} thickness={sl.thickness} nm" + f" enabled={sl.enabled}") + except Exception: + traceback.print_exc() + + print("\nsmoke probe DONE") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_buildup.py b/tests/test_buildup.py new file mode 100644 index 0000000..6d0faf7 --- /dev/null +++ b/tests/test_buildup.py @@ -0,0 +1,106 @@ +"""Solder/mask-opening buildup tests. Uniform-coverage and split-coverage +strips are 1D-exact, validating the harmonic-mean face weights and the +parallel-sheet conductance model to solver precision.""" +import numpy as np +import pytest + +from fill_resistance import raster, solver +from fill_resistance.geometry import (Polygon, SurfaceBuildup, load_problem, + save_problem) +from tests.util import NM, make_problem, ring_mm, sigma_s, strip_problem + +RHO_CU = 1.68e-8 +RHO_SN = 1.32e-7 + + +def _with_buildup(p, polys_mm, layer="F.Cu", solder_um=50.0, extra_um=0.0): + p.buildups = [SurfaceBuildup( + layer_name=layer, + polygons=[Polygon(outline=ring_mm(pts)) for pts in polys_mm])] + p.solder_thickness_nm = int(solder_um * 1000) + p.solder_rho_ohm_m = RHO_SN + p.extra_cu_nm = int(extra_um * 1000) + return p + + +def _solve(problem, h_mm): + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + return solver.run_solve(problem, stack, e1, e2, 1.0, + contact_model="equipotential"), stack + + +def _sigma_buildup(solder_um=50.0, extra_um=0.0): + return solder_um * 1e-6 / RHO_SN + extra_um * 1e-6 / RHO_CU + + +def test_full_coverage_exact(): + """Buildup over the whole strip: uniform parallel sheet, exact.""" + plain = strip_problem(length=50, width=10, e_len=5) + covered = _with_buildup(strip_problem(length=50, width=10, e_len=5), + [[(0, 0), (50, 0), (50, 10), (0, 10)]]) + r0, _ = _solve(plain, 0.5) + r1, stack = _solve(covered, 0.5) + sig = sigma_s() + _sigma_buildup() + assert r1.R_ohm == pytest.approx(81 / 20 / sig, rel=1e-9) + assert r1.R_ohm < r0.R_ohm + assert stack.buildup is not None and stack.buildup.any() + + +def test_half_coverage_series_exact(): + """Buildup over the right half: plain faces + one harmonic-mean + interface face + buildup faces in series, exact.""" + p = _with_buildup(strip_problem(length=50, width=10, e_len=5), + [[(25, 0), (50, 0), (50, 10), (25, 10)]]) + res, _ = _solve(p, 0.5) + s1 = sigma_s() + s2 = s1 + _sigma_buildup() + r_row = 40 / s1 + (s1 + s2) / (2 * s1 * s2) + 40 / s2 + assert res.R_ohm == pytest.approx(r_row / 20, rel=1e-9) + + +def test_buildup_only_over_hole_is_inert(): + """Solder wets copper only: an opening over a hole changes nothing.""" + holed = [([(0, 0), (50, 0), (50, 10), (0, 10)], + [[(20, 2), (30, 2), (30, 8), (20, 8)]])] + plain = make_problem(holed, rect1_mm=(0, 0, 5, 10), + rect2_mm=(45, 0, 50, 10)) + masked = _with_buildup( + make_problem(holed, rect1_mm=(0, 0, 5, 10), + rect2_mm=(45, 0, 50, 10)), + [[(21, 3), (29, 3), (29, 7), (21, 7)]]) # strictly inside the hole + r0, _ = _solve(plain, 0.5) + r1, _ = _solve(masked, 0.5) + assert r1.R_ohm == pytest.approx(r0.R_ohm, rel=1e-12) + + +def test_extra_copper_helps_more_than_solder(): + base = strip_problem(length=50, width=10, e_len=5) + solder_only = _with_buildup(strip_problem(length=50, width=10, e_len=5), + [[(0, 0), (50, 0), (50, 10), (0, 10)]]) + with_cu = _with_buildup(strip_problem(length=50, width=10, e_len=5), + [[(0, 0), (50, 0), (50, 10), (0, 10)]], + extra_um=70.0) + r0, _ = _solve(base, 0.5) + r1, _ = _solve(solder_only, 0.5) + r2, _ = _solve(with_cu, 0.5) + assert r2.R_ohm < r1.R_ohm < r0.R_ohm + # 50 um SAC solder ~ 6.4 um Cu: expect a modest (<15%) improvement + assert r1.R_ohm > 0.85 * r0.R_ohm + # +70 um Cu roughly halves R (2x thickness + solder) + assert r2.R_ohm < 0.55 * r0.R_ohm + + +def test_json_v4_roundtrip(tmp_path): + p = _with_buildup(strip_problem(), [[(0, 0), (50, 0), (50, 10), (0, 10)]], + extra_um=35.0) + f = tmp_path / "d.json" + save_problem(p, f) + q = load_problem(f) + assert len(q.buildups) == 1 and q.buildups[0].layer_name == "F.Cu" + assert q.solder_thickness_nm == 50_000 + assert q.extra_cu_nm == 35_000 + assert q.solder_rho_ohm_m == pytest.approx(RHO_SN) + r_p, _ = _solve(p, 0.5) + r_q, _ = _solve(q, 0.5) + assert r_q.R_ohm == pytest.approx(r_p.R_ohm, rel=1e-12) diff --git a/tests/test_contacts.py b/tests/test_contacts.py new file mode 100644 index 0000000..ea7193e --- /dev/null +++ b/tests/test_contacts.py @@ -0,0 +1,172 @@ +"""Contact-model (uniform vs equipotential) and multi-part terminal tests.""" +import numpy as np +import pytest + +from fill_resistance import raster, solver +from fill_resistance.errors import ElectrodeError +from fill_resistance.geometry import Electrode +from tests.util import NM, make_problem, rect_mm, sigma_s, strip_problem + + +def _solve(problem, h_mm, model, i_test=1.0): + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + parts1, parts2 = raster.electrode_partition(stack, problem) + return solver.run_solve(problem, stack, e1, e2, i_test, + contact_model=model, + parts1=parts1, parts2=parts2), stack + + +def _uniform_1d_reference(m_cols, n1, n2, sig, rows): + """Independent 1D reference: uniform injection over the first n1 + columns, extraction over the last n2, unit total current. R from + mean potentials, per row conductance sig, `rows` parallel rows.""" + inj = np.zeros(m_cols) + inj[:n1] += 1.0 / n1 + inj[m_cols - n2:] -= 1.0 / n2 + face_current = np.cumsum(inj)[:-1] # current through face k,k+1 + v = np.zeros(m_cols) + v[1:] = -np.cumsum(face_current) / sig # per single row of cells + v_plus = v[:n1].mean() + v_minus = v[m_cols - n2:].mean() + return (v_plus - v_minus) / 1.0 / rows # rows in parallel + + +def test_uniform_strip_exact_1d(): + """Full-width contacts on a uniform strip: rows are identical 1D + chains; compare with an independent 1D computation, exact.""" + p = strip_problem(length=50, width=10, e_len=5) + res, _ = _solve(p, 0.5, "uniform") + sig = sigma_s() + r_ref = _uniform_1d_reference(m_cols=100, n1=10, n2=10, sig=sig, rows=20) + assert res.R_ohm == pytest.approx(r_ref, rel=1e-9) + assert res.contact_model == "uniform" + assert res.power_balance_rel < 1e-9 # P = b^T V = I^2 R identity + + +def test_uniform_higher_than_equipotential(): + p = strip_problem(length=50, width=10, e_len=5) + r_uni, _ = _solve(p, 0.5, "uniform") + r_equ, _ = _solve(p, 0.5, "equipotential") + assert r_uni.R_ohm > r_equ.R_ohm + + +def test_uniform_current_density_ramps_inside_contact(): + """Inside the V+ contact, |J| must ramp: ~0 at the outer edge, + ~full sheet current at the inner (leading) edge; the equipotential + model shows ~0 throughout the contact interior.""" + p = strip_problem(length=50, width=10, e_len=5) + res_u, stack = _solve(p, 0.5, "uniform") + ny, nx = stack.shape2d + row = ny // 2 + # contact columns are the first 10 copper columns (margin = 2) + j_outer = res_u.Jmag[0, row, 2] # first contact column + j_inner = res_u.Jmag[0, row, 11] # last contact column + j_free = res_u.Jmag[0, row, nx // 2] # mid strip = I/(W t) + assert j_inner > 0.8 * j_free # ramped up to ~full + assert j_outer < 0.2 * j_free # near zero at outer edge + assert j_inner > 5 * max(j_outer, 1e-30) + + res_e, _ = _solve(p, 0.5, "equipotential") + j_center_e = res_e.Jmag[0, row, 6] # deep inside Dirichlet region + assert j_center_e < 0.05 * j_free + + +def test_multipart_terminal_equals_single_rect(): + """V+ split into two half-height rectangles == one full rectangle, + for both contact models (exact).""" + whole = strip_problem(length=50, width=10, e_len=5) + split = strip_problem(length=50, width=10, e_len=5) + split.electrodes1 = [ + Electrode(rect=rect_mm((0, 0, 5, 5))), + Electrode(rect=rect_mm((0, 5, 5, 10))), + ] + for model in ("uniform", "equipotential"): + r_whole, _ = _solve(whole, 0.5, model) + r_split, _ = _solve(split, 0.5, model) + assert r_split.R_ohm == pytest.approx(r_whole.R_ohm, rel=1e-9), model + + +def test_multipart_asymmetric_parts(): + """Two separated V+ parts feeding one V-: sane R, balance holds.""" + p = make_problem([([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 0, 2, 3), rect2_mm=(45, 0, 50, 10)) + p.electrodes1 = [ + Electrode(rect=rect_mm((0, 0, 2, 3)), label="top lug"), + Electrode(rect=rect_mm((0, 7, 2, 10)), label="bottom lug"), + ] + single, _ = _solve( + make_problem([([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 0, 2, 3), rect2_mm=(45, 0, 50, 10)), + 0.25, "uniform") + multi, _ = _solve(p, 0.25, "uniform") + assert multi.R_ohm < single.R_ohm # more contact area helps + assert multi.power_balance_rel < 1e-9 + + +def test_part_off_copper_raises_with_label(): + p = strip_problem(length=50, width=10, e_len=5) + p.electrodes1 = [ + Electrode(rect=rect_mm((0, 0, 5, 10))), + Electrode(rect=rect_mm((100, 100, 105, 105)), label="stray part"), + ] + stack = raster.rasterize_stack(p, 0.5 * NM) + with pytest.raises(ElectrodeError, match="stray part"): + raster.electrode_masks(stack, p) + + +def test_injection_area_currents_equipotential_flux(): + """Two V+ lugs at different distances: the nearer one carries more; + the flux split sums exactly to the test current.""" + p = make_problem([([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 0, 2, 10), rect2_mm=(48, 0, 50, 10)) + p.electrodes1 = [ + Electrode(rect=rect_mm((0, 4, 2, 6)), label="far lug"), + Electrode(rect=rect_mm((10, 4, 12, 6)), label="near lug"), + ] + res, _ = _solve(p, 0.25, "equipotential", i_test=10.0) + pc = dict(res.part_currents1) + assert pc["near lug"] > pc["far lug"] + assert pc["near lug"] + pc["far lug"] == pytest.approx(10.0, rel=1e-9) + # V- side: single part carries everything + assert res.part_currents2[0][1] == pytest.approx(10.0, rel=1e-9) + + +def test_injection_area_currents_uniform_area_share(): + """Uniform model: each injection area carries exactly its cell share.""" + p = make_problem([([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 0, 2, 10), rect2_mm=(48, 0, 50, 10)) + p.electrodes1 = [ + Electrode(rect=rect_mm((0, 0, 2, 6)), label="big"), # 2x6 mm + Electrode(rect=rect_mm((0, 6, 2, 9)), label="small"), # 2x3 mm + ] + res, _ = _solve(p, 0.25, "uniform", i_test=9.0) + pc = dict(res.part_currents1) + # cell counts: 8x24 = 192 and 8x12 = 96 at h=0.25 -> shares 2/3, 1/3 + assert pc["big"] == pytest.approx(6.0, rel=1e-12) + assert pc["small"] == pytest.approx(3.0, rel=1e-12) + + +def test_injection_area_partition_first_wins(): + """Overlapping parts: shared cells attributed to the first part, so + the shares still sum to the terminal current.""" + p = make_problem([([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 0, 2, 10), rect2_mm=(48, 0, 50, 10)) + p.electrodes1 = [ + Electrode(rect=rect_mm((0, 0, 2, 6)), label="first"), + Electrode(rect=rect_mm((0, 4, 2, 10)), label="second"), # overlaps + ] + res, _ = _solve(p, 0.25, "uniform", i_test=1.0) + total = sum(a for _, a in res.part_currents1) + assert total == pytest.approx(1.0, rel=1e-12) + + +def test_touching_ok_uniform_error_equipotential(): + p = make_problem([([(0, 0), (10, 0), (10, 10), (0, 10)], [])], + rect1_mm=(0, 0, 5, 10), rect2_mm=(5, 0, 10, 10)) + res, _ = _solve(p, 0.5, "uniform") # touching is fine here + assert res.R_ohm > 0 + p2 = make_problem([([(0, 0), (10, 0), (10, 10), (0, 10)], [])], + rect1_mm=(0, 0, 5, 10), rect2_mm=(5, 0, 10, 10)) + with pytest.raises(ElectrodeError, match="touch"): + _solve(p2, 0.5, "equipotential") \ No newline at end of file diff --git a/tests/test_geometry.py b/tests/test_geometry.py new file mode 100644 index 0000000..78bdee8 --- /dev/null +++ b/tests/test_geometry.py @@ -0,0 +1,88 @@ +import json +import math + +import numpy as np + +from fill_resistance.geometry import (arc_points, linearize_ring, + load_problem, problem_from_json, + save_problem) +from tests.util import make_multilayer, strip_problem + + +def test_arc_points_quarter_circle(): + r = 10_000_000 # 10 mm in nm + start = (r, 0) + mid = (int(r / math.sqrt(2)), int(r / math.sqrt(2))) + end = (0, r) + tol = 50_000 # 50 um + pts = arc_points(start, mid, end, tol) + assert len(pts) >= 3 + radii = np.hypot(pts[:, 0].astype(float), pts[:, 1].astype(float)) + assert np.allclose(radii, r, rtol=1e-6) + allpts = np.vstack([pts, [end]]).astype(float) + for a, b in zip(allpts[:-1], allpts[1:]): + half_chord = np.hypot(*(b - a)) / 2 + sagitta = r - math.sqrt(max(r**2 - half_chord**2, 0.0)) + assert sagitta <= tol * 1.01 + + +def test_arc_points_collinear_degrades_to_segment(): + pts = arc_points((0, 0), (5_000_000, 0), (10_000_000, 0), 1000) + assert len(pts) == 1 + assert tuple(pts[0]) == (0, 0) + + +def test_linearize_ring_mixed_nodes_and_closure(): + nodes = [ + ("pt", (0, 0)), + ("pt", (10, 0)), + ("arc", ((10, 0), (17, 7), (10, 14))), + ("pt", (0, 14)), + ("pt", (0, 0)), + ] + ring = linearize_ring(nodes, tol_nm=1) + assert (ring[0] != ring[-1]).any() + assert len(ring) > 4 + + +def test_problem_json_roundtrip_v2(tmp_path): + p = make_multilayer( + [[([(0, 0), (10, 0), (10, 1), (0, 1)], [])], + [([(0, 0), (10, 0), (10, 1), (0, 1)], [])]], + rect1_mm=(0, 0, 1, 1), rect2_mm=(9, 0, 10, 1), + contact1="L0", contact2="L1", vias_mm=[(5.5, 0.5)]) + f = tmp_path / "dump.json" + save_problem(p, f) + q = load_problem(f) + assert q.layer_names == ["L0", "L1"] + assert q.electrodes1[0].contact == "L0" + assert q.electrodes2[0].contact == "L1" + assert len(q.vias) == 1 and q.vias[0].drill_nm == p.vias[0].drill_nm + assert q.plating_nm == p.plating_nm + assert np.array_equal(q.layers[0].polygons[0].outline, + p.layers[0].polygons[0].outline) + assert abs(q.sigma_s(0) - p.sigma_s(0)) < 1e-12 * p.sigma_s(0) + + +def test_v1_schema_still_loads(): + p = strip_problem() + v1 = { + "schema_version": 1, + "board_path": "old", + "layer_name": "F.Cu", + "net_name": "GND", + "thickness_nm": 70000, + "thickness_source": "stackup", + "rho_ohm_m": 1.68e-8, + "rect1": {"x0": 0, "y0": 0, "x1": 1000, "y1": 1000, + "layer_name": "User.1"}, + "rect2": {"x0": 5000, "y0": 0, "x1": 6000, "y1": 1000, + "layer_name": "User.1"}, + "polygons": [{"outline": p.layers[0].polygons[0].outline.tolist(), + "holes": []}], + } + q = problem_from_json(json.loads(json.dumps(v1))) + assert q.layer_names == ["F.Cu"] + assert q.vias == [] + assert q.electrodes1[0].contact == "all" + assert q.layers[0].thickness_nm == 70000 diff --git a/tests/test_multilayer.py b/tests/test_multilayer.py new file mode 100644 index 0000000..3835eba --- /dev/null +++ b/tests/test_multilayer.py @@ -0,0 +1,148 @@ +"""Multi-layer / via solver tests. The 1-cell-wide strip cases are pure +series chains, so the discrete solution is exact and validates the via +barrel model, layer coupling, and flux integration to solver precision.""" +import math + +import numpy as np +import pytest + +from fill_resistance import raster, solver +from fill_resistance.errors import ConnectivityError, ElectrodeError +from tests.util import NM, make_multilayer, sigma_s + + +def _solve(problem, h_mm, i_test=1.0): + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + return solver.run_solve(problem, stack, e1, e2, i_test, + contact_model="equipotential"), stack + + +STRIP = [(0, 0), (10, 0), (10, 1), (0, 1)] # 10 x 1 mm; h=1 -> 1 row + + +def _r_via(length_mm, drill_mm=0.3, plating_um=18.0, rho=1.68e-8): + area = math.pi * (drill_mm * 1e-3) * (plating_um * 1e-6) + return rho * (length_mm * 1e-3) / area + + +def test_two_identical_layers_parallel(): + """Both electrodes contact both layers, no vias needed: two identical + independent sheets in parallel -> exactly half the single-layer R.""" + one = make_multilayer([[(STRIP, [])]], (0, 0, 1, 1), (9, 0, 10, 1)) + two = make_multilayer([[(STRIP, [])], [(STRIP, [])]], + (0, 0, 1, 1), (9, 0, 10, 1)) + r1, _ = _solve(one, 1.0) + r2, _ = _solve(two, 1.0) + assert r2.R_ohm == pytest.approx(r1.R_ohm / 2, rel=1e-9) + + +def test_via_chain_1d_exact(): + """e1 on L0 left end, e2 on L1 right end, one through-via at x=5.5: + R = faces_L0/sigma + R_via + faces_L1/sigma, exact.""" + p = 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) + res, stack = _solve(p, 1.0) + # cells at x centers 0.5..9.5 -> cols 0..9 (plus margins); electrode + # cells: col of 0.5 (e1), col of 9.5 (e2); via cell: col of 5.5 + sig = sigma_s() + faces_l0 = 5 # cols 0->5: faces between 0|1 .. 4|5 + faces_l1 = 4 # cols 5->9 + r_exact = (faces_l0 + faces_l1) / sig + _r_via(1.0) + assert res.R_ohm == pytest.approx(r_exact, rel=1e-9) + assert res.mismatch_rel < 1e-10 + assert len(res.via_reports) == 1 + # the single via carries the full test current + assert res.via_reports[0].current_a == pytest.approx(1.0, rel=1e-9) + assert res.via_reports[0].power_w == pytest.approx(_r_via(1.0), rel=1e-9) + + +def test_parallel_vias_halve_barrel_resistance(): + base = dict(rect1_mm=(0, 0, 1, 1), rect2_mm=(9, 0, 10, 1), + contact1="L0", contact2="L1", gap_mm=1.0) + one = make_multilayer([[(STRIP, [])], [(STRIP, [])]], + vias_mm=[(5.5, 0.5)], **base) + two = make_multilayer([[(STRIP, [])], [(STRIP, [])]], + vias_mm=[(5.5, 0.5), (5.5, 0.5)], **base) + r_one, _ = _solve(one, 1.0) + r_two, _ = _solve(two, 1.0) + assert (r_one.R_ohm - r_two.R_ohm) == pytest.approx(_r_via(1.0) / 2, + rel=1e-9) + + +def test_antipad_bridging(): + """3 layers; the middle layer has an antipad hole at the via cell, so + the barrel bridges L0 -> L2 directly with DOUBLE the length.""" + mid_with_hole = [(STRIP, [[(5, 0), (6, 0), (6, 1), (5, 1)]])] + p = make_multilayer( + [[(STRIP, [])], mid_with_hole, [(STRIP, [])]], + rect1_mm=(0, 0, 1, 1), rect2_mm=(9, 0, 10, 1), + contact1="L0", contact2="L2", + vias_mm=[(5.5, 0.5)], gap_mm=1.0) + res, _ = _solve(p, 1.0) + sig = sigma_s() + r_exact = (5 + 4) / sig + _r_via(2.0) # barrel length 2 mm + assert res.R_ohm == pytest.approx(r_exact, rel=1e-9) + + +def test_via_short_between_electrodes_raises(): + """Both electrodes over the SAME cell on different layers with a via + there = direct short, no free copper -> error.""" + p = make_multilayer( + [[(STRIP, [])], [(STRIP, [])]], + rect1_mm=(5, 0, 6, 1), rect2_mm=(5, 0, 6, 1), + contact1="L0", contact2="L1", + vias_mm=[(5.5, 0.5)], gap_mm=1.0) + stack = raster.rasterize_stack(p, 1.0 * NM) + e1, e2 = raster.electrode_masks(stack, p) + with pytest.raises(ElectrodeError, match="directly connected"): + solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + + +def test_layers_without_via_disconnected(): + p = make_multilayer( + [[(STRIP, [])], [(STRIP, [])]], + rect1_mm=(0, 0, 1, 1), rect2_mm=(9, 0, 10, 1), + contact1="L0", contact2="L1", gap_mm=1.0) # no vias + stack = raster.rasterize_stack(p, 1.0 * NM) + e1, e2 = raster.electrode_masks(stack, p) + with pytest.raises(ConnectivityError, match="not connected"): + solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + + +def test_power_split_layers_and_vias(): + """Power accounting: layer + via powers sum to I^2 R.""" + p = 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) + res, _ = _solve(p, 1.0, i_test=10.0) + assert res.power_balance_rel < 1e-9 + assert res.P_vias == pytest.approx(100 * _r_via(1.0), rel=1e-9) + sig = sigma_s() + assert res.P_layers[0] == pytest.approx(100 * 5 / sig, rel=1e-9) + assert res.P_layers[1] == pytest.approx(100 * 4 / sig, rel=1e-9) + + +def test_pad_polygon_electrode(): + """A polygon-shaped electrode (pad) restricted to one layer.""" + from fill_resistance.geometry import Electrode, Polygon + from tests.util import rect_mm, ring_mm + p = make_multilayer([[(STRIP, [])]], (0, 0, 1, 1), (9, 0, 10, 1)) + # replace terminal 1 with a small polygon pad covering the same cell + p.electrodes1 = [Electrode( + rect=rect_mm((0.2, 0.2, 0.8, 0.8)), + contact="all", + polygons=[Polygon(outline=ring_mm( + [(0.2, 0.2), (0.8, 0.2), (0.8, 0.8), (0.2, 0.8)]))], + label="pad TP1.1", + )] + res, _ = _solve(p, 1.0) + sig = sigma_s() + assert res.R_ohm == pytest.approx(9 / sig, rel=1e-9) # cols 0..9 diff --git a/tests/test_raster.py b/tests/test_raster.py new file mode 100644 index 0000000..7737a19 --- /dev/null +++ b/tests/test_raster.py @@ -0,0 +1,102 @@ +import numpy as np +import pytest + +from fill_resistance import config, raster, solver +from fill_resistance.errors import (ConnectivityError, ElectrodeError, + GridSizeError) +from tests.util import NM, make_problem, strip_problem + + +def _stack(problem, h_mm): + return raster.rasterize_stack(problem, h_mm * NM) + + +def test_exact_cell_count_square_with_hole(): + # 10x10 mm square, centered 4x4 mm hole, h=1 mm: cell centers at + # half-integers, no boundary ambiguity -> exactly 100 - 16 cells + p = make_problem( + [([(0, 0), (10, 0), (10, 10), (0, 10)], + [[(3, 3), (7, 3), (7, 7), (3, 7)]])], + rect1_mm=(0, 0, 1, 10), rect2_mm=(9, 0, 10, 10)) + stack = _stack(p, 1.0) + assert int(stack.masks[0].sum()) == 100 - 16 + + +def test_margin_cells_are_empty(): + p = strip_problem() + stack = _stack(p, 0.5) + m = stack.masks[0] + assert not m[0, :].any() and not m[-1, :].any() + assert not m[:, 0].any() and not m[:, -1].any() + + +def test_electrode_masks_and_counts(): + p = strip_problem(length=50, width=10, e_len=5) + stack = _stack(p, 0.5) + e1, e2 = raster.electrode_masks(stack, p) + # electrodes: 5 mm / 0.5 mm = 10 columns x 20 rows on 1 layer + assert int(e1.sum()) == 200 and int(e2.sum()) == 200 + + +def test_electrode_off_copper_raises(): + p = make_problem([([(0, 0), (10, 0), (10, 10), (0, 10)], [])], + rect1_mm=(20, 20, 25, 25), rect2_mm=(8, 0, 10, 10)) + stack = _stack(p, 0.5) + with pytest.raises(ElectrodeError, match="does not overlap"): + raster.electrode_masks(stack, p) + + +def test_touching_electrodes_raise_equipotential(): + p = make_problem([([(0, 0), (10, 0), (10, 10), (0, 10)], [])], + rect1_mm=(0, 0, 5, 10), rect2_mm=(5, 0, 10, 10)) + stack = _stack(p, 0.5) + e1, e2 = raster.electrode_masks(stack, p) # touch is fine at mask level + with pytest.raises(ElectrodeError, match="touch"): + solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + + +def test_disconnected_regions_raise(): + p = make_problem( + [([(0, 0), (10, 0), (10, 10), (0, 10)], []), + ([(20, 0), (30, 0), (30, 10), (20, 10)], [])], + rect1_mm=(0, 0, 2, 10), rect2_mm=(28, 0, 30, 10)) + stack = _stack(p, 0.5) + e1, e2 = raster.electrode_masks(stack, p) + with pytest.raises(ConnectivityError, match="not connected"): + solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + + +def test_islands_dropped(): + # island square not touching the main strip disappears from the mask + p = make_problem( + [([(0, 0), (50, 0), (50, 10), (0, 10)], []), + ([(20, 20), (30, 20), (30, 30), (20, 30)], [])], + rect1_mm=(0, 0, 5, 10), rect2_mm=(45, 0, 50, 10)) + stack = _stack(p, 0.5) + e1, e2 = raster.electrode_masks(stack, p) + before = int(stack.masks.sum()) + solver.run_solve(p, stack, e1, e2, 1.0, contact_model="equipotential") + after = int(stack.masks.sum()) + assert after < before + assert after == 100 * 20 # only the strip remains + + +def test_hard_max_cells_guard(monkeypatch): + monkeypatch.setattr(config, "CELL_UM_OVERRIDE", 1.0) # 1 um cells + p = strip_problem() + with pytest.raises(GridSizeError, match="M cells"): + raster.choose_cell_size(p.copper_bbox(), len(p.layers)) + + +def test_auto_cell_size_hits_target(): + # large plane: unclamped regime, cell count tracks TARGET_CELLS + p = strip_problem(length=200, width=100, e_len=5) + h = raster.choose_cell_size(p.copper_bbox(), 1) + ncells = (200.0 * NM / h) * (100.0 * NM / h) + assert 0.5 * config.TARGET_CELLS < ncells < 2.0 * config.TARGET_CELLS + # small board: MIN_CELL_UM clamp kicks in instead + q = strip_problem() # 50x10 mm + hq = raster.choose_cell_size(q.copper_bbox(), 1) + assert hq == pytest.approx(config.MIN_CELL_UM * 1000) diff --git a/tests/test_skin.py b/tests/test_skin.py new file mode 100644 index 0000000..d5d6394 --- /dev/null +++ b/tests/test_skin.py @@ -0,0 +1,103 @@ +"""Skin-effect model tests. The single-layer AC solve scales ALL in-plane +conductances identically, so R_AC = R_DC * resistance_factor EXACTLY - +which turns the analytic foil formula into an end-to-end exact test.""" +import math + +import numpy as np +import pytest + +from fill_resistance import raster, skin, solver +from tests.util import NM, make_multilayer, strip_problem + +RHO = 1.68e-8 + + +def _solve(problem, h_mm, i_test=1.0, freq=0.0): + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + return solver.run_solve(problem, stack, e1, e2, i_test, freq, + contact_model="equipotential"), stack + + +def test_skin_depth_value(): + # copper @ 1 MHz: ~65-66 um + assert skin.skin_depth_m(1e6, RHO) * 1e6 == pytest.approx(65.2, rel=0.01) + + +def test_sheet_resistance_dc_limit(): + t = 70e-6 + assert skin.sheet_resistance_ac(t, 0.0, RHO) == pytest.approx(RHO / t) + # low frequency: within 0.1% of DC + assert skin.sheet_resistance_ac(t, 100.0, RHO) == pytest.approx( + RHO / t, rel=1e-3) + + +def test_sheet_resistance_high_f_limits(): + t = 70e-6 + f = 1e9 # delta << t + delta = skin.skin_depth_m(f, RHO) + assert skin.sheet_resistance_ac(t, f, RHO, sides=1) == pytest.approx( + RHO / delta, rel=0.01) + assert skin.sheet_resistance_ac(t, f, RHO, sides=2) == pytest.approx( + RHO / (2 * delta), rel=0.01) + + +def test_resistance_factor_monotonic(): + t = 70e-6 + factors = [skin.resistance_factor(t, f, RHO) + for f in (0, 1e4, 1e5, 1e6, 1e7, 1e8)] + assert all(b >= a - 1e-12 for a, b in zip(factors, factors[1:])) + assert factors[0] == 1.0 + + +def test_parse_frequency(): + assert skin.parse_frequency("") == 0.0 + assert skin.parse_frequency("0") == 0.0 + assert skin.parse_frequency("142k") == 142_000.0 + assert skin.parse_frequency("1.5M") == 1_500_000.0 + assert skin.parse_frequency("2meg") == 2_000_000.0 + assert skin.parse_frequency("100000") == 100_000.0 + assert skin.parse_frequency("100 kHz") == 100_000.0 + assert skin.parse_frequency("junk") == 0.0 + + +def test_single_layer_ac_scales_exactly(): + """Uniform conductance scaling leaves the field shape unchanged: + R_AC = R_DC * factor to solver precision.""" + p = strip_problem(length=50, width=10, e_len=5) + f = 2e6 # delta=46um < t=70um + r_dc, _ = _solve(p, 0.5) + r_ac, _ = _solve(p, 0.5, freq=f) + factor = skin.resistance_factor(70e-6, f, RHO, sides=1) + assert factor > 1.2 # real crowding at 2 MHz + assert r_ac.R_ohm == pytest.approx(r_dc.R_ohm * factor, rel=1e-9) + assert r_ac.rs_ratios[0] == pytest.approx(factor, rel=1e-12) + assert r_ac.skin_depth_um == pytest.approx( + skin.skin_depth_m(f, RHO) * 1e6, rel=1e-12) + + +def test_via_chain_ac_exact(): + """1D chain: layers scale by the foil factor, the barrel by the + plating-wall factor - exact composition.""" + STRIP = [(0, 0), (10, 0), (10, 1), (0, 1)] + p = 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) + f = 2e6 + sig_ac = 1.0 / skin.sheet_resistance_ac(70e-6, f, RHO, sides=1) + via_factor = skin.resistance_factor(18e-6, f, RHO, sides=2) + r_via_dc = RHO * 1e-3 / (math.pi * 0.3e-3 * 18e-6) + r_exact = (5 + 4) / sig_ac + r_via_dc * via_factor + res, _ = _solve(p, 1.0, freq=f) + assert res.R_ohm == pytest.approx(r_exact, rel=1e-9) + + +def test_dc_default_unchanged(): + """freq omitted -> identical to the pre-skin behavior.""" + p = strip_problem(length=50, width=10, e_len=5) + res, _ = _solve(p, 0.5) + assert res.freq_hz == 0.0 + assert res.skin_depth_um is None + assert all(r == 1.0 for r in res.rs_ratios) diff --git a/tests/test_solver.py b/tests/test_solver.py new file mode 100644 index 0000000..bf8a4db --- /dev/null +++ b/tests/test_solver.py @@ -0,0 +1,124 @@ +import numpy as np +import pytest + +from fill_resistance import config, raster, solver +from tests.util import NM, make_problem, strip_problem + + +def _solve(problem, h_mm, i_test=1.0): + stack = raster.rasterize_stack(problem, h_mm * NM) + e1, e2 = raster.electrode_masks(stack, problem) + return solver.run_solve(problem, stack, e1, e2, i_test, + contact_model="equipotential"), stack + + +def test_uniform_strip_exact_discrete(): + """Full-width electrodes on a uniform strip: every row is an identical + series chain, so the discrete solution is exact: + R = (n_free_columns + 1) / (n_rows * sigma_s).""" + p = strip_problem(length=50, width=10, e_len=5) + res, stack = _solve(p, 0.5) + n_free_cols = int(round((50 - 2 * 5) / 0.5)) # 80 + n_rows = int(round(10 / 0.5)) # 20 + r_exact = (n_free_cols + 1) / n_rows / p.sigma_s(0) + assert res.R_ohm == pytest.approx(r_exact, rel=1e-9) + assert res.mismatch_rel < 1e-10 + r_cont = p.rho_ohm_m * 0.0405 / (0.010 * 70e-6) + assert res.R_ohm == pytest.approx(r_cont, rel=1e-6) + + +def test_strip_R_independent_of_h(): + p = strip_problem(length=50, width=10, e_len=5) + for h in (1.0, 0.5, 0.25): + res, _ = _solve(p, h) + m = int(round(40 / h)) + rows = int(round(10 / h)) + assert res.R_ohm == pytest.approx((m + 1) / rows / p.sigma_s(0), + rel=1e-9) + + +def test_partial_electrode_constriction(): + full = strip_problem(length=50, width=10, e_len=5) + partial = make_problem( + [([(0, 0), (50, 0), (50, 10), (0, 10)], [])], + rect1_mm=(0, 4, 5, 6), + rect2_mm=(45, 4, 50, 6)) + r_full, _ = _solve(full, 0.25) + r_a, _ = _solve(partial, 0.25) + r_b, _ = _solve(partial, 0.125) + assert r_a.R_ohm > r_full.R_ohm * 1.05 + assert abs(r_a.R_ohm - r_b.R_ohm) < 0.01 * r_b.R_ohm + + +def test_l_shape_corner_squares(): + """Right-angle bend of equal-width arms: corner square counts as + ~0.559 squares (conformal-mapping result), within 5% at a fine grid.""" + w = 10.0 + outline = [(0, 0), (30, 0), (30, 30), (20, 30), (20, 10), (0, 10)] + p = make_problem([(outline, [])], + rect1_mm=(0, 0, 2, 10), + rect2_mm=(20, 28, 30, 30)) + res, _ = _solve(p, 0.125) + a_sq = (20 - 2) / w + b_sq = (28 - 10) / w + r_expect = (a_sq + b_sq + 0.559) / p.sigma_s(0) + assert res.R_ohm == pytest.approx(r_expect, rel=0.05) + + +def test_hole_increases_resistance_and_converges(): + solid = make_problem([([(0, 0), (40, 0), (40, 20), (0, 20)], [])], + rect1_mm=(0, 0, 2, 20), rect2_mm=(38, 0, 40, 20)) + holed = make_problem( + [([(0, 0), (40, 0), (40, 20), (0, 20)], + [[(15, 5), (25, 5), (25, 15), (15, 15)]])], + rect1_mm=(0, 0, 2, 20), rect2_mm=(38, 0, 40, 20)) + r_solid, _ = _solve(solid, 0.25) + r_a, _ = _solve(holed, 0.25) + r_b, _ = _solve(holed, 0.125) + assert r_a.R_ohm > r_solid.R_ohm * 1.1 + assert abs(r_a.R_ohm - r_b.R_ohm) < 0.01 * r_b.R_ohm + + +def test_cg_path_matches_direct(monkeypatch): + p = strip_problem(length=50, width=10, e_len=5) + r_direct, _ = _solve(p, 0.25) + monkeypatch.setattr(config, "SPSOLVE_MAX_UNKNOWNS", 0) + r_cg, _ = _solve(p, 0.25) + assert r_cg.solve_info.method == "cg+jacobi" + assert r_cg.R_ohm == pytest.approx(r_direct.R_ohm, rel=1e-6) + assert r_cg.mismatch_rel < 1e-5 + + +def test_current_density_and_potential_scale(): + """Uniform strip at 1 A: |J| in the free region equals I/(W t); the + potential span equals R * I.""" + p = strip_problem(length=50, width=10, e_len=5) + res, stack = _solve(p, 0.5) + j_expect = 1.0 / (0.010 * 70e-6) + ny, nx = stack.shape2d + assert res.Jmag[0, ny // 2, nx // 2] == pytest.approx(j_expect, rel=1e-6) + assert np.nanmax(res.V) == pytest.approx(res.R_ohm, rel=1e-9) + + +def test_test_current_scaling(): + """V and J scale linearly with I_test, power quadratically; R fixed.""" + p = strip_problem(length=50, width=10, e_len=5) + r1, _ = _solve(p, 0.5, i_test=1.0) + r10, _ = _solve(p, 0.5, i_test=10.0) + assert r10.R_ohm == pytest.approx(r1.R_ohm, rel=1e-12) + assert np.nanmax(r10.V) == pytest.approx(10 * np.nanmax(r1.V), rel=1e-9) + assert np.nanmax(r10.Jmag) == pytest.approx(10 * np.nanmax(r1.Jmag), + rel=1e-9) + assert r10.P_total == pytest.approx(100 * r1.P_total, rel=1e-9) + + +def test_power_identity(): + """Sum of edge powers equals I^2 R exactly for the direct solve.""" + p = make_problem( + [([(0, 0), (40, 0), (40, 20), (0, 20)], + [[(15, 5), (25, 5), (25, 15), (15, 15)]])], + rect1_mm=(0, 0, 2, 20), rect2_mm=(38, 0, 40, 20)) + res, _ = _solve(p, 0.25, i_test=40.0) + assert res.power_balance_rel < 1e-9 + assert res.P_total == pytest.approx(40.0 ** 2 * res.R_ohm, rel=1e-12) + assert res.P_vias == 0.0 diff --git a/tests/util.py b/tests/util.py new file mode 100644 index 0000000..be7dbcb --- /dev/null +++ b/tests/util.py @@ -0,0 +1,74 @@ +"""Synthetic Problem builders for tests. Dimensions in mm here, nm inside.""" +import numpy as np + +from fill_resistance.geometry import (Electrode, LayerFill, Polygon, Problem, + Rect, ViaLink) + +NM = 1_000_000 # nm per mm + + +def ring_mm(points_mm) -> np.ndarray: + return (np.asarray(points_mm, dtype=float) * NM).astype(np.int64) + + +def _polys(polygons_mm) -> list[Polygon]: + return [ + Polygon(outline=ring_mm(outline), holes=[ring_mm(h) for h in holes]) + for outline, holes in polygons_mm + ] + + +def rect_mm(r, layer="User.1") -> Rect: + return Rect.normalized(int(r[0] * NM), int(r[1] * NM), + int(r[2] * NM), int(r[3] * NM), layer) + + +def make_multilayer(layers_mm, rect1_mm, rect2_mm, contact1="all", + contact2="all", vias_mm=(), t_um=70.0, gap_mm=1.0, + drill_mm=0.3, plating_um=18.0, rho=1.68e-8) -> Problem: + """layers_mm: one list of (outline_pts, [holes]) per layer; layer i is + named 'L{i}' at z = i * gap_mm. vias_mm: (x, y) through-barrels.""" + nlayers = len(layers_mm) + return Problem( + board_path="synthetic", + net_name="TEST", + rho_ohm_m=rho, + plating_nm=int(plating_um * 1000), + layers=[ + LayerFill(layer_name=f"L{i}", thickness_nm=int(t_um * 1000), + z_nm=int(i * gap_mm * NM), polygons=_polys(polys_mm)) + for i, polys_mm in enumerate(layers_mm) + ], + vias=[ + ViaLink(x=int(x * NM), y=int(y * NM), + drill_nm=int(drill_mm * NM), + z_top_nm=-1, z_bot_nm=int((nlayers - 1) * gap_mm * NM) + 1) + for x, y in vias_mm + ], + electrodes1=[Electrode(rect=rect_mm(rect1_mm), contact=contact1)], + electrodes2=[Electrode(rect=rect_mm(rect2_mm), contact=contact2)], + thickness_source="override", + ) + + +def make_problem(polygons_mm, rect1_mm, rect2_mm, t_um=70.0, + rho=1.68e-8) -> Problem: + """Single-layer problem (the v1 test surface).""" + p = make_multilayer([polygons_mm], rect1_mm, rect2_mm, t_um=t_um, rho=rho) + p.layers[0].layer_name = "F.Cu" + return p + + +def strip_problem(length=50.0, width=10.0, e_len=5.0, t_um=70.0): + """Uniform strip with full-width electrodes at both ends.""" + outline = [(0, 0), (length, 0), (length, width), (0, width)] + return make_problem( + [(outline, [])], + rect1_mm=(0, 0, e_len, width), + rect2_mm=(length - e_len, 0, length, width), + t_um=t_um, + ) + + +def sigma_s(t_um=70.0, rho=1.68e-8) -> float: + return t_um * 1e-6 / rho