diff --git a/README.md b/README.md index 505cde6..715db6f 100644 --- a/README.md +++ b/README.md @@ -114,9 +114,10 @@ SWIG API. Requires KiCad **10.0.1+**. 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 (~2 M cells total - across layers by default). Direct sparse solve up to 2.5 M unknowns, - Jacobi-CG above. Discretization error typically ≲ 2 % at defaults — - halve the cell size and compare to judge convergence. + across layers by default). Direct sparse solve up to 500 k unknowns, + AMG-preconditioned CG (pyamg) above — Jacobi-CG if pyamg is missing. + Discretization error typically ≲ 2 % at defaults — halve the cell size + and compare to judge convergence. ## Offline / development @@ -133,7 +134,7 @@ Linux/macOS use `.venv/bin/python`): ```powershell uv venv --python 3.11 .venv -uv pip install --python .venv\Scripts\python.exe kicad-python numpy scipy matplotlib pytest +uv pip install --python .venv\Scripts\python.exe kicad-python numpy scipy pyamg matplotlib pytest .venv\Scripts\python.exe -m pytest tests -q # incl. exact analytic cases .venv\Scripts\python.exe tools\api_probe.py # IPC API probe vs live KiCad .venv\Scripts\python.exe -m fill_resistance.board_io dump.json [NET] # extract only diff --git a/fill_resistance/config.py b/fill_resistance/config.py index e52b602..f3bced3 100644 --- a/fill_resistance/config.py +++ b/fill_resistance/config.py @@ -5,10 +5,9 @@ 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. +# 150 um cells down to 50 um. 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 @@ -51,10 +50,13 @@ CONTACT_MODEL = "uniform" # "uniform": conductor pressed on top injects # (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 +SPSOLVE_MAX_UNKNOWNS = 500_000 # above this, AMG-preconditioned CG (pyamg; + # Jacobi-CG if pyamg is missing). Direct is + # exact and fastest for small grids; measured + # at 1.4M unknowns: spsolve 13 s / ~3 GB, + # AMG-CG 6 s at a fraction of the memory +AMG_TOL = 1e-10 # relative residual of the AMG-CG solve +CG_TOL = 1e-8 # Jacobi-CG fallback (no pyamg) CG_MAXITER = 50_000 # CG iterations are cheap; large grids need many # --- Geometry --- diff --git a/fill_resistance/raster.py b/fill_resistance/raster.py index d10770d..19672ad 100644 --- a/fill_resistance/raster.py +++ b/fill_resistance/raster.py @@ -18,6 +18,7 @@ from dataclasses import dataclass import numpy as np from matplotlib.path import Path as MplPath +from PIL import Image, ImageDraw from scipy import ndimage from . import config @@ -106,7 +107,10 @@ def choose_cell_size(bbox_nm: tuple[int, int, int, int], nlayers: int) -> float: 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).""" + working only within the ring's bbox. Hybrid rasterizer: PIL scanline + fill for the bulk (fast, O(vertices + cells)), then the cells within + a ~2 px band around the ring edge are re-tested exactly against the + polygon, so the result is identical to a pure center-in-polygon pass.""" ny, nx = stack.shape2d h = stack.h_nm j0 = max(0, int((ring[:, 0].min() - stack.x0_nm) / h) - 1) @@ -115,13 +119,35 @@ def _paint_ring(stack: RasterStack, ring: np.ndarray, value: bool, 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) + w, ht = j1 - j0, i1 - i0 + + # cell (i, j) center <-> pixel (j - j0, i - i0) + px = (ring[:, 0] - stack.x0_nm) / h - 0.5 - j0 + py = (ring[:, 1] - stack.y0_nm) / h - 0.5 - i0 + pts = list(zip(px.tolist(), py.tolist())) + + inside = np.zeros((ht, w), dtype=bool) + band = np.ones((ht, w), dtype=bool) + if len(pts) >= 3: + fill_img = Image.new("1", (w, ht), 0) + ImageDraw.Draw(fill_img).polygon(pts, fill=1) + inside = np.array(fill_img, dtype=bool) + band_img = Image.new("1", (w, ht), 0) + ImageDraw.Draw(band_img).line(pts + pts[:1], fill=1, width=5, + joint="curve") + band = np.array(band_img, dtype=bool) + + bi, bj = np.nonzero(band) + if len(bi): + xs = stack.x0_nm + (bj + j0 + 0.5) * h + ys = stack.y0_nm + (bi + i0 + 0.5) * h + # 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[bi, bj] = MplPath(verts, closed=True).contains_points( + np.column_stack([xs, ys])) + sub = target[i0:i1, j0:j1] sub[inside] = value diff --git a/fill_resistance/solver.py b/fill_resistance/solver.py index cbb1974..0d6d850 100644 --- a/fill_resistance/solver.py +++ b/fill_resistance/solver.py @@ -293,10 +293,40 @@ def solve_system(A: sparse.csr_matrix, b: np.ndarray) -> tuple[np.ndarray, Solve if n <= config.SPSOLVE_MAX_UNKNOWNS: x = sla.spsolve(A.tocsc(), b) return x, SolveInfo(method="spsolve", n_unknowns=n) + try: + return _solve_amg(A, b) + except ImportError: + print("note: pyamg not installed - falling back to Jacobi-CG " + "(much slower on large grids)") + return _solve_cg_jacobi(A, b) + +def _solve_amg(A: sparse.csr_matrix, b: np.ndarray) -> tuple[np.ndarray, SolveInfo]: + """CG preconditioned with smoothed-aggregation AMG: near-linear + scaling on these 2D Laplacians and a fraction of spsolve's memory.""" + import pyamg + + n = A.shape[0] + ml = pyamg.smoothed_aggregation_solver(A.tocsr(), max_coarse=500) + residuals: list[float] = [] + x = ml.solve(b, tol=config.AMG_TOL, maxiter=300, accel="cg", + residuals=residuals) + res = float(np.linalg.norm(b - A @ x) / max(np.linalg.norm(b), 1e-300)) + if not np.isfinite(res) or res > 1e-6: + raise SolverError( + f"AMG-CG did not converge (residual {res:.2e}). Try a " + f"different grid size, or force the direct solver by raising " + f"SPSOLVE_MAX_UNKNOWNS in config.py." + ) + return x, SolveInfo(method="amg+cg", n_unknowns=n, + iterations=max(len(residuals) - 1, 0), residual=res) + + +def _solve_cg_jacobi(A: sparse.csr_matrix, b: np.ndarray) -> tuple[np.ndarray, SolveInfo]: # 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. + n = A.shape[0] d = A.diagonal() M = sla.LinearOperator((n, n), lambda v: v / d) iters = 0 diff --git a/fill_resistance/standalone.py b/fill_resistance/standalone.py index 0813a56..4999e31 100644 --- a/fill_resistance/standalone.py +++ b/fill_resistance/standalone.py @@ -42,7 +42,8 @@ def main(argv=None) -> int: 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") + help="use the iterative solver (AMG-CG, or Jacobi-CG " + "without pyamg) regardless of problem size") args = ap.parse_args(argv) if args.cell_um is not None: diff --git a/metadata.json b/metadata.json index 419d45e..160b64d 100644 --- a/metadata.json +++ b/metadata.json @@ -2,7 +2,7 @@ "$schema": "https://go.kicad.org/pcm/schemas/v2", "name": "Fill Resistance", "description": "DC/AC resistance of copper zone fills between two contacts, single- or multi-layer with via coupling, with current and power density maps.", - "description_full": "Computes the DC or AC resistance of copper zone fills between two contacts (marker rectangles on User.1/User.2 and/or selected pads), single- or multi-layer: the chosen net's fills are solved as coupled finite-difference sheets linked by the net's via and through-hole-pad barrels.\n\nShows per-layer rasterized maps, potential, current density and power density, reports per-via currents (via ampacity) and total dissipation at a selectable test current. At a user-set frequency the exact 1D foil/barrel skin-effect correction is applied (AC results are a rigorous lower bound). PNGs, a text summary and a re-solvable geometry dump are saved per run.\n\nNote: the first load builds the plugin's Python environment (numpy, scipy, matplotlib, PySide6) and can take several minutes.", + "description_full": "Computes the DC or AC resistance of copper zone fills between two contacts (marker rectangles on User.1/User.2 and/or selected pads), single- or multi-layer: the chosen net's fills are solved as coupled finite-difference sheets linked by the net's via and through-hole-pad barrels.\n\nShows per-layer rasterized maps, potential, current density and power density, reports per-via currents (via ampacity) and total dissipation at a selectable test current. At a user-set frequency the exact 1D foil/barrel skin-effect correction is applied (AC results are a rigorous lower bound). PNGs, a text summary and a re-solvable geometry dump are saved per run.\n\nNote: the first load builds the plugin's Python environment (numpy, scipy, pyamg, matplotlib, PySide6) and can take several minutes.", "identifier": "th.co.b4l.fill-resistance", "type": "plugin", "author": { diff --git a/requirements.txt b/requirements.txt index 919a38e..da3eb5b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ kicad-python>=0.7.0 numpy scipy +pyamg matplotlib PySide6 diff --git a/tests/test_raster.py b/tests/test_raster.py index d305104..c156983 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -22,6 +22,35 @@ def test_exact_cell_count_square_with_hole(): assert int(stack.masks[0].sum()) == 100 - 16 +def test_hybrid_raster_matches_exact_point_test(): + """The PIL-fill + exact-edge-band rasterizer must be cell-for-cell + identical to a pure center-in-polygon pass, including awkward + fractional offsets, concave lobes and a hole.""" + from matplotlib.path import Path as MplPath + ang = np.linspace(0, 2 * np.pi, 257, endpoint=False) + r = 7.3 + 1.7 * np.sin(5 * ang) + 0.9 * np.cos(9 * ang + 0.4) + blob = np.stack([20.05 + r * np.cos(ang), 20.13 + r * np.sin(ang)], + axis=1) + hole = np.stack([20.4 + 2.1 * np.cos(ang), 19.8 + 2.2 * np.sin(ang)], + axis=1) + p = make_problem([(blob.tolist(), [hole.tolist()])], + rect1_mm=(14, 19, 16, 21), rect2_mm=(24, 19, 26, 21)) + stack = _stack(p, 0.25) + + ny, nx = stack.shape2d + xg, yg = stack.cell_centers(0, ny, 0, nx) + pts = np.column_stack([xg.ravel(), yg.ravel()]) + + def exact(ring): + verts = np.vstack([ring, ring[:1]]) + return MplPath(verts, closed=True).contains_points(pts).reshape( + ny, nx) + + poly = p.layers[0].polygons[0] + ref = exact(poly.outline) & ~exact(poly.holes[0]) + assert np.array_equal(stack.masks[0], ref) + + def test_margin_cells_are_empty(): p = strip_problem() stack = _stack(p, 0.5) diff --git a/tests/test_solver.py b/tests/test_solver.py index b929bf5..9898869 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -79,10 +79,21 @@ def test_hole_increases_resistance_and_converges(): assert abs(r_a.R_ohm - r_b.R_ohm) < 0.01 * r_b.R_ohm -def test_cg_path_matches_direct(monkeypatch): +def test_iterative_paths_match_direct(monkeypatch): + """AMG-CG (default iterative) and Jacobi-CG (pyamg-missing fallback) + both reproduce the direct solve.""" p = strip_problem(length=50, width=10, e_len=5) r_direct, _ = _solve(p, 0.25) monkeypatch.setattr(config, "SPSOLVE_MAX_UNKNOWNS", 0) + r_amg, _ = _solve(p, 0.25) + assert r_amg.solve_info.method == "amg+cg" + assert r_amg.R_ohm == pytest.approx(r_direct.R_ohm, rel=1e-6) + assert r_amg.mismatch_rel < 1e-5 + + def no_pyamg(A, b): + raise ImportError("pyamg unavailable") + + monkeypatch.setattr(solver, "_solve_amg", no_pyamg) 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)