From 14cdf4c7ee8d9bacb7e325bd202b1e659c3190d8 Mon Sep 17 00:00:00 2001 From: janik Date: Wed, 15 Jul 2026 17:05:31 +0700 Subject: [PATCH] Add balanced quadtree grid engine (adaptive cells, phase 1) fill_resistance/quadtree.py decomposes a layer's fine copper mask into 2:1-balanced power-of-two leaves: boundary and keep-fine cells stay at the fine size, interiors coarsen with their Chebyshev distance to the nearest feature (guard factor, default 4), and an explicit enforcement pass splits any leaf more than twice an edge-adjacent neighbor. Face conductances use the series-half-cell rule, which reduces to the production harmonic mean for equal sizes and EXACTLY to sigma in the uniform limit - verified edge-for-edge against solver.build_edges and to rel 1e-12 in R against run_solve, so the exact-value test suite stays authoritative for this engine. Measured (feature-dense 120x120 plate, h=50um, production AMG solver): uniform 5.58M unknowns ~35s; adaptive guard=4 823k / ~8s at -1.1%; guard=8 1.74M / ~12s at -0.47%. tools/adaptive_proto.py now benchmarks the engine itself. Not yet wired into the pipeline: phase 2 ports electrodes, barrels, 1D chains, buildup and field output onto leaves. Co-Authored-By: Claude Fable 5 --- fill_resistance/quadtree.py | 171 ++++++++++++++++++++++++++++++++++++ tests/test_quadtree.py | 141 +++++++++++++++++++++++++++++ tools/adaptive_proto.py | 169 ++++++++++------------------------- 3 files changed, 357 insertions(+), 124 deletions(-) create mode 100644 fill_resistance/quadtree.py create mode 100644 tests/test_quadtree.py diff --git a/fill_resistance/quadtree.py b/fill_resistance/quadtree.py new file mode 100644 index 0000000..1126ec5 --- /dev/null +++ b/fill_resistance/quadtree.py @@ -0,0 +1,171 @@ +"""Adaptive quadtree grid engine (phase 1 of the adaptive-cell work). + +Decomposes a layer's fine copper mask into a 2:1-BALANCED set of square +leaves (power-of-two sizes in fine-cell units): cells at copper +boundaries and in keep_fine regions stay at the fine size, interiors +coarsen with their Chebyshev distance to the nearest feature (guard +factor), and an explicit enforcement pass splits any leaf more than +twice the size of an edge-adjacent neighbor. + +The uniform limit (max_block=1) reproduces the fine grid EXACTLY: one +leaf per copper cell and face conductance sigma - the production +solver's grid is the special case, which keeps the exact-value test +suite authoritative for this engine. + +Face conductance between edge-adjacent leaves a, b sharing w fine-cell +faces is the series-half-cell expression + + g = w / (size_a / (2 sigma_a) + size_b / (2 sigma_b)) + +which reduces to the harmonic mean 2 sigma_a sigma_b / (sigma_a + +sigma_b) for equal sizes (the production buildup/mouth face rule) and +to sigma in the uniform limit. + +Phase 2 (not here) wires this into the pipeline: electrodes, barrels, +1D chains, buildup and field output still run on the uniform grid. +""" +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy import ndimage + + +@dataclass +class LeafGrid: + """Square leaves over one layer, in fine-cell units.""" + y0: np.ndarray # int32, aligned: y0 % size == 0 + x0: np.ndarray + size: np.ndarray # int32, power of two + id_grid: np.ndarray # (ny, nx) int32; -1 = not copper + + @property + def n(self) -> int: + return len(self.size) + + +def build_leaves(mask: np.ndarray, keep_fine: np.ndarray | None = None, + max_block: int = 32, guard: int = 4) -> LeafGrid: + """Balanced leaf decomposition of a boolean copper mask. keep_fine + marks cells that must stay at the fine size (electrodes, and later + via mouths / buildup edges). guard scales how much clearance a block + of size s needs (Chebyshev distance >= guard * s).""" + ny, nx = mask.shape + if keep_fine is None: + keep_fine = np.zeros_like(mask) + coarsenable = mask & ~keep_fine + + # Chebyshev distance to the nearest non-coarsenable cell (array + # border padded as background so edges never look like interior) + pad = np.pad(coarsenable, 1) + D = ndimage.distance_transform_cdt(pad, metric="chessboard")[1:-1, 1:-1] + + S = np.zeros((ny, nx), dtype=np.int32) + S[mask] = 1 + s = 2 + while s <= max_block: + S[D >= guard * s] = s + s *= 2 + + grid = _emit(S, mask, max_block) + for _ in range(32): + if _split_unbalanced(grid, S): + grid = _emit(S, mask, max_block) + else: + return grid + raise RuntimeError("quadtree balance did not converge") + + +def _emit(S: np.ndarray, mask: np.ndarray, max_block: int) -> LeafGrid: + """Greedy top-down emission: an aligned s-block becomes a leaf where + every cell allows size s and nothing larger claimed it.""" + ny, nx = mask.shape + py, px = (-ny) % max_block, (-nx) % max_block + Sp = np.pad(S, ((0, py), (0, px))) + NY, NX = Sp.shape + id_grid = np.full((NY, NX), -1, dtype=np.int32) + covered = np.zeros((NY, NX), dtype=bool) + y0s, x0s, sizes = [], [], [] + nid = 0 + s = max_block + while s >= 2: + min_S = Sp.reshape(NY // s, s, NX // s, s).min(axis=(1, 3)) + free = ~covered.reshape(NY // s, s, NX // s, s).any(axis=(1, 3)) + cand = (min_S >= s) & free + k = int(cand.sum()) + if k: + lvl = np.full(cand.shape, -1, dtype=np.int32) + lvl[cand] = nid + np.arange(k, dtype=np.int32) + up = np.repeat(np.repeat(lvl, s, axis=0), s, axis=1) + sel = up >= 0 + id_grid[sel] = up[sel] + covered |= sel + ii, jj = np.nonzero(cand) + y0s.append(ii.astype(np.int32) * s) + x0s.append(jj.astype(np.int32) * s) + sizes.append(np.full(k, s, dtype=np.int32)) + nid += k + s //= 2 + fi, fj = np.nonzero(np.pad(mask, ((0, py), (0, px))) & ~covered) + id_grid[fi, fj] = nid + np.arange(len(fi), dtype=np.int32) + y0s.append(fi.astype(np.int32)) + x0s.append(fj.astype(np.int32)) + sizes.append(np.ones(len(fi), dtype=np.int32)) + return LeafGrid( + y0=np.concatenate(y0s) if y0s else np.zeros(0, np.int32), + x0=np.concatenate(x0s) if x0s else np.zeros(0, np.int32), + size=np.concatenate(sizes) if sizes else np.zeros(0, np.int32), + id_grid=id_grid[:ny, :nx], + ) + + +def _split_unbalanced(grid: LeafGrid, S: np.ndarray) -> bool: + """Cap S over any leaf more than 2x an edge-adjacent neighbor, so the + next emission splits it. Returns True if anything was capped.""" + ia, ib, _ = leaf_faces(grid) + sa, sb = grid.size[ia], grid.size[ib] + big = np.unique(np.concatenate([ia[sa > 2 * sb], ib[sb > 2 * sa]])) + for lid in big: + y, x, s = int(grid.y0[lid]), int(grid.x0[lid]), int(grid.size[lid]) + np.minimum(S[y:y + s, x:x + s], s // 2, out=S[y:y + s, x:x + s]) + return len(big) > 0 + + +def leaf_faces(grid: LeafGrid) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """(ia, ib, w) for every edge-adjacent leaf pair, each pair once; + w = shared face length in fine-cell units.""" + n = grid.n + if n == 0: + z = np.zeros(0, dtype=np.int64) + return z, z, z + out = [] + for sl_a, sl_b in ((np.s_[:, :-1], np.s_[:, 1:]), + (np.s_[:-1, :], np.s_[1:, :])): + a = grid.id_grid[sl_a].ravel().astype(np.int64) + b = grid.id_grid[sl_b].ravel().astype(np.int64) + ok = (a >= 0) & (b >= 0) & (a != b) + key, counts = np.unique(a[ok] * n + b[ok], return_counts=True) + out.append((key // n, key % n, counts)) + ia = np.concatenate([o[0] for o in out]) + ib = np.concatenate([o[1] for o in out]) + w = np.concatenate([o[2] for o in out]) + return ia, ib, w + + +def leaf_edges(grid: LeafGrid, sigma) -> tuple[np.ndarray, np.ndarray, + np.ndarray]: + """Face conductances [S]: sigma is a scalar or per-leaf array of + sheet conductance. Uniform limit -> exactly sigma per face.""" + ia, ib, w = leaf_faces(grid) + sig = np.broadcast_to(np.asarray(sigma, dtype=float), (grid.n,)) + g = w / (grid.size[ia] / (2.0 * sig[ia]) + + grid.size[ib] / (2.0 * sig[ib])) + return ia, ib, g + + +def balanced(grid: LeafGrid) -> bool: + """2:1 balance invariant: edge-adjacent leaves differ <= 2x in size.""" + ia, ib, _ = leaf_faces(grid) + sa, sb = grid.size[ia], grid.size[ib] + return bool((np.maximum(sa, sb) <= 2 * np.minimum(sa, sb)).all()) diff --git a/tests/test_quadtree.py b/tests/test_quadtree.py new file mode 100644 index 0000000..ba4408b --- /dev/null +++ b/tests/test_quadtree.py @@ -0,0 +1,141 @@ +"""Quadtree grid engine tests (phase 1): exact uniform limit against the +production graph and solver, partition/alignment/balance invariants, and +adaptive-vs-fine R agreement.""" +import numpy as np +import pytest +from scipy import ndimage + +from fill_resistance import quadtree, raster, solver +from tests.util import NM, make_problem, strip_problem + + +def _plate_with_holes(n=5, size_mm=40.0): + holes = [] + pitch = size_mm / n + for i in range(n): + for j in range(n): + x, y = pitch * (i + 0.4), pitch * (j + 0.4) + holes.append([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)]) + outline = [(0, 0), (size_mm, 0), (size_mm, size_mm), (0, size_mm)] + return make_problem([(outline, holes)], + rect1_mm=(0, 15, 2, 25), + rect2_mm=(size_mm - 2, 15, size_mm, 25)) + + +def _solve_on_leaves(problem, stack, e1, e2, grid): + """Equipotential mini-solve on the leaf graph, using the production + assembly and linear solver.""" + ia, ib, g = quadtree.leaf_edges(grid, problem.sigma_s(0)) + state = np.ones(grid.n, dtype=np.uint8) + for e, code in ((e1[0], 2), (e2[0], 3)): + ids = grid.id_grid[e] + state[ids[ids >= 0]] = code + edges = solver.Edges(a=ia, b=ib, w=g, + via_index=np.full(len(ia), -1, dtype=np.int32)) + A, rhs, _ = solver._assemble(state, edges, None) + x, _ = solver.solve_system(A, rhs) + V = np.zeros(grid.n) + V[state == 2] = 1.0 + V[state == 1] = x + Ie = g * (V[ia] - V[ib]) + sa, sb = state[ia], state[ib] + I1 = float(Ie[sa == 2].sum() - Ie[sb == 2].sum()) + I2 = float(Ie[sb == 3].sum() - Ie[sa == 3].sum()) + return 1.0 / (0.5 * (I1 + I2)) + + +def test_uniform_limit_graph_identical(): + """max_block=1: one leaf per cell, and the edge list matches the + production in-plane graph exactly (same pairs, conductance sigma).""" + p = _plate_with_holes() + stack = raster.rasterize_stack(p, 0.5 * NM) + grid = quadtree.build_leaves(stack.masks[0], max_block=1) + assert int(stack.masks.sum()) == grid.n + assert (grid.size == 1).all() + + ny, nx = stack.shape2d + flat_of_leaf = grid.y0.astype(np.int64) * nx + grid.x0 + ia, ib, g = quadtree.leaf_edges(grid, p.sigma_s(0)) + ours = np.sort(np.stack([ + np.minimum(flat_of_leaf[ia], flat_of_leaf[ib]), + np.maximum(flat_of_leaf[ia], flat_of_leaf[ib])], axis=1), axis=0) + + edges = solver.build_edges(stack, p, [p.sigma_s(0)]) + ref = np.sort(np.stack([np.minimum(edges.a, edges.b), + np.maximum(edges.a, edges.b)], axis=1), axis=0) + assert ours.shape == ref.shape + assert np.array_equal(np.sort(ours.view("i8,i8"), order=["f0", "f1"], + axis=0), + np.sort(ref.view("i8,i8"), order=["f0", "f1"], + axis=0)) + assert np.allclose(g, p.sigma_s(0), rtol=0, atol=0) + + +def test_uniform_limit_R_matches_production(): + p = strip_problem(length=50, width=10, e_len=5) + stack = raster.rasterize_stack(p, 0.25 * NM) + e1, e2 = raster.electrode_masks(stack, p) + ref = solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + stack2 = raster.rasterize_stack(p, 0.25 * NM) + e1b, e2b = raster.electrode_masks(stack2, p) + grid = quadtree.build_leaves(stack2.masks[0], max_block=1) + R = _solve_on_leaves(p, stack2, e1b, e2b, grid) + assert R == pytest.approx(ref.R_ohm, rel=1e-12) + + +def test_partition_alignment_and_balance(): + p = _plate_with_holes() + stack = raster.rasterize_stack(p, 0.1 * NM) + mask = stack.masks[0] + grid = quadtree.build_leaves(mask, max_block=32) + + # exact partition of the copper + assert int((grid.size.astype(np.int64) ** 2).sum()) == int(mask.sum()) + assert (grid.id_grid >= 0).sum() == int(mask.sum()) + assert not (grid.id_grid[~mask] >= 0).any() + counts = np.bincount(grid.id_grid[grid.id_grid >= 0], minlength=grid.n) + assert np.array_equal(counts, grid.size.astype(np.int64) ** 2) + + # power-of-two sizes, aligned to their own size + assert np.array_equal(grid.size & (grid.size - 1), + np.zeros_like(grid.size)) + assert (grid.y0 % grid.size == 0).all() + assert (grid.x0 % grid.size == 0).all() + + # 2:1 balance and real coarsening (max size is geometry-limited by + # the guard distance, not by max_block, on this feature-dense plate) + assert quadtree.balanced(grid) + assert grid.n < 0.5 * int(mask.sum()) + assert grid.size.max() >= 4 + + +def test_boundary_and_keep_fine_stay_fine(): + p = _plate_with_holes() + stack = raster.rasterize_stack(p, 0.1 * NM) + mask = stack.masks[0] + keep = np.zeros_like(mask) + keep[50:60, 50:60] = True + grid = quadtree.build_leaves(mask, keep_fine=keep, max_block=32) + + boundary = mask & ndimage.binary_dilation(~mask) + assert (grid.size[grid.id_grid[boundary]] == 1).all() + assert (grid.size[grid.id_grid[keep & mask]] == 1).all() + + +def test_adaptive_R_close_to_fine(): + """Adaptive leaves reproduce the fine-uniform R within 1% on the + holey plate (features everywhere - the adversarial case).""" + p = _plate_with_holes() + stack = raster.rasterize_stack(p, 0.1 * NM) + e1, e2 = raster.electrode_masks(stack, p) + ref = solver.run_solve(p, stack, e1, e2, 1.0, + contact_model="equipotential") + + stack2 = raster.rasterize_stack(p, 0.1 * NM) + e1b, e2b = raster.electrode_masks(stack2, p) + grid = quadtree.build_leaves(stack2.masks[0], + keep_fine=(e1b[0] | e2b[0])) + R = _solve_on_leaves(p, stack2, e1b, e2b, grid) + assert grid.n < 0.4 * ref.solve_info.n_unknowns + assert R == pytest.approx(ref.R_ohm, rel=0.01) diff --git a/tools/adaptive_proto.py b/tools/adaptive_proto.py index bc0b26e..6ca5f32 100644 --- a/tools/adaptive_proto.py +++ b/tools/adaptive_proto.py @@ -1,27 +1,25 @@ -"""Prototype: quadtree-adaptive grid for the fill-resistance solver. +"""Benchmark for the adaptive quadtree grid engine +(fill_resistance/quadtree.py, phase 1 of the adaptive-cell work). -Coarsens the existing fine raster bottom-up (power-of-two blocks that are -fully copper away from electrodes, erosion-graded per level), builds the -leaf graph fully vectorized via an id-grid (face conductance -g = sigma * overlap / mean-size, which reduces EXACTLY to the production -sigma in the uniform limit), and reuses the production assembly + AMG -solver. Run from the repo root: .venv/Scripts/python tools/adaptive_proto.py + .venv/Scripts/python tools/adaptive_proto.py -Measured 2026-07-15 (120x120 mm plate, 400 holes, h = 50 um; adversarial: -features everywhere, so geometric refinement has no smooth interior): +Solves a feature-dense 120x120 mm plate (400 holes) at h = 50 um on the +uniform grid and on the balanced quadtree (engine defaults: guard=4, +max_block=32), using the production assembly + solver for both. - uniform R = 0.508504 mOhm 5.58M unknowns 27 s (reference) - max_block=4 R = 0.506368 mOhm 474k unknowns 3 s -0.42% 12x - max_block=8 R = 0.503386 mOhm 253k unknowns 1 s -1.0% 22x - max_block=16 R = 0.497873 mOhm 213k unknowns 1 s -2.1% 26x +Measured 2026-07-15 on this machine: + + uniform R = 0.508504 mOhm 5.58M unknowns ~30-40 s + adaptive guard=4 R = 0.502793 mOhm 823k unknowns ~8 s -1.1% + adaptive guard=8 R = 0.506102 mOhm 1.74M unknowns ~12 s -0.47% uniform-limit check (max_block=1): rel diff 0.00e+00 vs production. -Coarsening biases R low (coarse cells overestimate conductance where the -field curves); a production version needs true 2:1 balancing + a guard -band, and optionally one residual-driven refine pass, to push the -max_block=4 accuracy to larger blocks. On big-pour boards (smooth -interiors) the unknown ratios are far higher than on this geometry. +This is the ADVERSARIAL case (features at 6 mm pitch everywhere, no +smooth interior); big-pour boards coarsen far more aggressively. The +residual bias is the first-order coarse-fine interface flux - phase 4 +(gradient-corrected fluxes / solution-adaptive refinement) is the +lever if tighter accuracy per leaf is needed. """ import sys import time @@ -31,127 +29,49 @@ import numpy as np sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) -from fill_resistance import raster, solver +from fill_resistance import quadtree, raster, solver from tests.util import NM, make_problem, strip_problem -MAX_BLOCK = 32 # coarsest leaf = 32 x 32 fine cells - -def build_leaves(mask, e1, e2, max_block=MAX_BLOCK): - """Greedy top-down coarsening. Returns (y0, x0, size) per leaf plus - an id_grid at fine resolution (-1 = empty).""" - ny, nx = mask.shape - pad_y = (-ny) % max_block - pad_x = (-nx) % max_block - m = np.pad(mask, ((0, pad_y), (0, pad_x))) - coars = m & ~np.pad(e1 | e2, ((0, pad_y), (0, pad_x))) - NY, NX = m.shape - - from scipy import ndimage - id_grid = np.full((NY, NX), -1, dtype=np.int64) - covered = np.zeros((NY, NX), dtype=bool) - y0s, x0s, sizes = [], [], [] - nid = 0 - levels = [] - red = coars.copy() - s = 1 - while s < max_block: - red = red.reshape(red.shape[0] // 2, 2, red.shape[1] // 2, 2 - ).all(axis=(1, 3)) - s *= 2 - # grading: a size-s block must sit in an all-copper 3x3 block - # neighborhood at its own level, so leaf sizes step down smoothly - # toward boundaries (guard band + approximate 2:1 balance) - graded = ndimage.binary_erosion(red, np.ones((3, 3), dtype=bool)) - levels.append((s, graded)) - for s, allc in reversed(levels): - cov_k = covered.reshape(NY // s, s, NX // s, s).any(axis=(1, 3)) - cand = allc & ~cov_k - ii, jj = np.nonzero(cand) - for i_, j_ in zip(ii, jj): - id_grid[i_ * s:(i_ + 1) * s, j_ * s:(j_ + 1) * s] = nid - y0s.append(i_ * s) - x0s.append(j_ * s) - sizes.append(s) - nid += 1 - covered |= np.repeat(np.repeat(cand, s, axis=0), s, axis=1) - fi, fj = np.nonzero(m & ~covered) - n_fine = len(fi) - id_grid[fi, fj] = nid + np.arange(n_fine) - y0s.extend(fi.tolist()) - x0s.extend(fj.tolist()) - sizes.extend([1] * n_fine) - return (np.array(y0s), np.array(x0s), np.array(sizes), - id_grid[:ny, :nx]) - - -def leaf_edges(id_grid, sizes, sigma): - """All leaf-leaf face conductances, vectorized: count shared fine - faces per leaf pair (= overlap length w), g = sigma * w / mean(sa, sb). - Uniform limit: w = 1, sizes 1 -> g = sigma (identical to production).""" - aa, bb, ww = [], [], [] - n = len(sizes) - for sl_a, sl_b in ((np.s_[:, :-1], np.s_[:, 1:]), - (np.s_[:-1, :], np.s_[1:, :])): - a = id_grid[sl_a].ravel() - b = id_grid[sl_b].ravel() - ok = (a >= 0) & (b >= 0) & (a != b) - key = a[ok] * n + b[ok] - uniq, counts = np.unique(key, return_counts=True) - ia = uniq // n - ib = uniq % n - g = sigma * counts / (0.5 * (sizes[ia] + sizes[ib])) - aa.append(ia) - bb.append(ib) - ww.append(g) - return (np.concatenate(aa), np.concatenate(bb), np.concatenate(ww)) - - -def solve_adaptive(problem, h_mm, max_block=MAX_BLOCK): +def solve_adaptive(problem, h_mm, **kw): stack = raster.rasterize_stack(problem, h_mm * NM) e1, e2 = raster.electrode_masks(stack, problem) t0 = time.perf_counter() - y0, x0, sizes, id_grid = build_leaves(stack.masks[0], e1[0], e2[0], - max_block) - a, b, w = leaf_edges(id_grid, sizes, problem.sigma_s(0)) + grid = quadtree.build_leaves(stack.masks[0], + keep_fine=(e1[0] | e2[0]), **kw) + ia, ib, g = quadtree.leaf_edges(grid, problem.sigma_s(0)) t_build = time.perf_counter() - t0 - n = len(sizes) - state = np.ones(n, dtype=np.uint8) - fine_ids = id_grid[e1[0]] - state[fine_ids[fine_ids >= 0]] = 2 - fine_ids = id_grid[e2[0]] - state[fine_ids[fine_ids >= 0]] = 3 - - edges = solver.Edges(a=a, b=b, w=w, - via_index=np.full(len(a), -1, dtype=np.int32)) + state = np.ones(grid.n, dtype=np.uint8) + for e, code in ((e1[0], 2), (e2[0], 3)): + ids = grid.id_grid[e] + state[ids[ids >= 0]] = code + edges = solver.Edges(a=ia, b=ib, w=g, + via_index=np.full(len(ia), -1, dtype=np.int32)) t0 = time.perf_counter() - A, rhs, idx = solver._assemble(state, edges, None) + A, rhs, _ = solver._assemble(state, edges, None) x, info = solver.solve_system(A, rhs) t_solve = time.perf_counter() - t0 - V = np.zeros(n) + V = np.zeros(grid.n) V[state == 2] = 1.0 V[state == 1] = x - Ie = w * (V[a] - V[b]) - sa, sb = state[a], state[b] + Ie = g * (V[ia] - V[ib]) + sa, sb = state[ia], state[ib] I1 = float(Ie[sa == 2].sum() - Ie[sb == 2].sum()) I2 = float(Ie[sb == 3].sum() - Ie[sa == 3].sum()) - R = 1.0 / (0.5 * (I1 + I2)) - mismatch = abs(I1 - I2) / max(abs(I1), abs(I2)) - return R, n, info, t_build, t_solve, mismatch + return 1.0 / (0.5 * (I1 + I2)), grid, info, t_build, t_solve -# --- 1) uniform-limit correctness: max_block=1 must equal production --- +# --- 1) uniform-limit correctness --- p = strip_problem(length=50, width=10, e_len=5) stack = raster.rasterize_stack(p, 0.25 * NM) e1, e2 = raster.electrode_masks(stack, p) res = solver.run_solve(p, stack, e1, e2, 1.0, contact_model="equipotential") -R_u, n_u, *_ = solve_adaptive(p, 0.25, max_block=1) -print(f"uniform-limit check: production R={res.R_ohm:.12g}, " - f"prototype R={R_u:.12g}, rel diff {abs(R_u / res.R_ohm - 1):.2e}") +R_u, *_ = solve_adaptive(p, 0.25, max_block=1) +print(f"uniform-limit check: rel diff {abs(R_u / res.R_ohm - 1):.2e}") -# --- 2) the payoff case: 120x120 plate, 400 holes, h = 50 um --- +# --- 2) feature-dense plate at h = 50 um --- holes = [] for i in range(20): for j in range(20): @@ -167,13 +87,14 @@ e1, e2 = raster.electrode_masks(stack, big) ref = solver.run_solve(big, stack, e1, e2, 1.0, contact_model="equipotential") t_ref = time.perf_counter() - t0 -print(f"\nuniform 50um : R = {ref.R_ohm * 1e3:.6g} mOhm, " - f"{ref.solve_info.n_unknowns} unknowns, {t_ref:.1f} s total " +print(f"uniform : R = {ref.R_ohm * 1e3:.6g} mOhm, " + f"{ref.solve_info.n_unknowns} unknowns, {t_ref:.1f} s " f"({ref.solve_info.method})") -for mb in (4, 8, 16, 32): - R_a, n_a, info, t_b, t_s, mm = solve_adaptive(big, 0.05, max_block=mb) - print(f"adaptive mb={mb:2d}: R = {R_a * 1e3:.6g} mOhm, " - f"{info.n_unknowns:8d} unknowns, build {t_b:.1f} s + " - f"solve {t_s:.1f} s, rel diff {abs(R_a / ref.R_ohm - 1):.2e}, " - f"ratio {ref.solve_info.n_unknowns / info.n_unknowns:.0f}x") +R_a, grid, info, t_b, t_s = solve_adaptive(big, 0.05) +print(f"adaptive: R = {R_a * 1e3:.6g} mOhm, {info.n_unknowns} unknowns, " + f"build {t_b:.1f} s + solve {t_s:.1f} s ({info.method})") +print(f"rel diff {abs(R_a / ref.R_ohm - 1):.2e}, " + f"{ref.solve_info.n_unknowns / info.n_unknowns:.1f}x fewer unknowns, " + f"max leaf {int(grid.size.max())} cells, balanced=" + f"{quadtree.balanced(grid)}")