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)}")