conduction.py

1import math2import sys3from fractions import Fraction4import time56import numpy as np7import scipy.sparse as sp8import scipy.sparse.linalg as spl910R3 = math.sqrt(3.0)11SIGMA = 1.0 / R31213# RENDER1415def odd_digits(n, level, k=1):16    i = np.arange(n ** level)17    out = []18    for _ in range(level):19        out.append((i % n) % 2 == 1)20        i = i // n21    out = np.array(out)22    if k > 1:23        out = np.repeat(out, k, axis=1)24    return out2526def render(n, level, k=1):27    o = odd_digits(n, level, k)28    void = np.zeros((o.shape[1], o.shape[1]), bool)29    for row in o:30        void |= np.outer(row, row)31    return ~void3233def obnosov(z):34    return math.sqrt((1 + 3 * z) / (3 + z))3536def fill(n):37    return n * n - ((n - 1) // 2) ** 23839# NETWORK4041def pair(a, b):42    s = a + b43    return np.where(s > 0, 2 * a * b / np.where(s > 0, s, 1), 0.0)4445def solve(n_nodes, rows, cols, w, diag, rhs):46    A = sp.coo_matrix((-w, (rows, cols)), shape=(n_nodes, n_nodes))47    A = (A + A.T + sp.diags(diag)).tocsc()48    return spl.spsolve(A, rhs, permc_spec="COLAMD")4950def network(a, wx, wy, left, right, right_value):51    nx, ny = a.shape52    live = a > 053    idx = -np.ones(a.shape, np.int64)54    idx[live] = np.arange(live.sum())55    n_nodes = int(live.sum())56    rows, cols, ws = [], [], []57    gx = pair(a[:-1, :], a[1:, :]) * wx[:-1, :]58    gy = pair(a[:, :-1], a[:, 1:]) * wy[:, :-1]59    for g, i0, i1 in ((gx, idx[:-1, :], idx[1:, :]), (gy, idx[:, :-1], idx[:, 1:])):60        m = g > 061        rows.append(i0[m])62        cols.append(i1[m])63        ws.append(g[m])64    rows, cols, ws = np.concatenate(rows), np.concatenate(cols), np.concatenate(ws)65    diag = np.bincount(rows, ws, n_nodes) + np.bincount(cols, ws, n_nodes)66    rhs = np.zeros(n_nodes)67    lm = live[0, :]68    diag[idx[0, lm]] += left[lm]69    rm = live[-1, :] & (right > 0)70    diag[idx[-1, rm]] += right[rm]71    rhs[idx[-1, rm]] += right[rm] * right_value72    u = solve(n_nodes, rows, cols, ws, diag, rhs)73    return float(np.sum(left[lm] * u[idx[0, lm]]))7475def conductance_full(a):76    a = np.asarray(a, float)77    one = np.ones(a.shape)78    return network(a, one, one, 2 * a[0, :], 2 * a[-1, :], 1.0)7980def conductance(a):81    a = np.asarray(a, float)82    s = a.shape[0]83    h = (s + 1) // 284    q = a[:h, :h].copy()85    wx = np.ones(q.shape)86    wy = np.ones(q.shape)87    rw = np.ones(h)88    if s % 2:89        wx[:, -1] = 0.590        rw[-1] = 0.591        center = q[-1, :] > 092        left = 2 * q[0, :] * rw93        g = pair(q[-2, :], q[-1, :]) * rw94        q2 = q[:-1, :]95        wx2, wy2 = wx[:-1, :], wy[:-1, :]96        right = np.where(center, g, 0.0)97        return 2 * network(q2, wx2, wy2, left, right, 0.5)98    left = 2 * q[0, :]99    right = 2 * q[-1, :]100    return 2 * network(q, wx, wy, left, right, 0.5)101102def rich(g1, g2, r=2.0, p=4.0 / 3.0):103    return g2 + (g2 - g1) / (r ** p - 1)104105def aitken(g1, g2, g3, r=2.0):106    e1, e2 = g2 - g1, g3 - g2107    p = math.log(e1 / e2) / math.log(r)108    return rich(g2, g3, r, p), p109110def level(n, lv, k):111    return conductance(render(n, lv, k).astype(float))112113# CELL114115def cell(m, z):116    a = np.ones((2 * m, 2 * m))117    a[m:, m:] = z118    return conductance_full(a)119120def cell3(m, rtol=1e-11):121    t = np.arange(2 * m) >= m122    up = t[:, None, None].astype(int) + t[None, :, None] + t[None, None, :]123    live = up <= 1124    s = 2 * m125    idx = -np.ones(live.shape, np.int64)126    n_nodes = int(live.sum())127    idx[live] = np.arange(n_nodes)128    rows, cols = [], []129    for ax in range(3):130        a = [slice(None)] * 3131        b = [slice(None)] * 3132        a[ax] = slice(0, -1)133        b[ax] = slice(1, None)134        both = live[tuple(a)] & live[tuple(b)]135        rows.append(idx[tuple(a)][both])136        cols.append(idx[tuple(b)][both])137    rows, cols = np.concatenate(rows), np.concatenate(cols)138    w = np.ones(len(rows))139    diag = np.bincount(rows, w, n_nodes) + np.bincount(cols, w, n_nodes)140    left = idx[0][live[0]]141    right = idx[-1][live[-1]]142    diag[left] += 2.0143    diag[right] += 2.0144    rhs = np.zeros(n_nodes)145    rhs[right] = 2.0146    A = sp.coo_matrix((-w, (rows, cols)), shape=(n_nodes, n_nodes))147    A = (A + A.T + sp.diags(diag)).tocsr()148    x0 = (np.indices(live.shape)[0][live] + 0.5) / s149    u, info = spl.cg(A, rhs, x0=x0, M=sp.diags(1.0 / diag), rtol=rtol, maxiter=20000)150    assert info == 0151    return 2.0 * float(np.sum(u[left])) / s152153def verb_cell():154    ms = (64, 128, 256, 512)155    g = [cell(m, 0.0) for m in ms]156    for i, m in enumerate(ms):157        line = f"insulating M={m} sigma={g[i]:.10f} err={g[i] - SIGMA:+.3e}"158        if i:159            line += f" order={math.log((g[i - 1] - SIGMA) / (g[i] - SIGMA)) / math.log(2):.4f} rich={rich(g[i - 1], g[i]):.10f}"160        print(line)161    e1, e2 = rich(g[1], g[2]), rich(g[2], g[3])162    bar = abs(e2 - e1)163    print(f"insulating sigma={e2:.10f} bar={bar:.1e} target={SIGMA:.10f} off={e2 - SIGMA:+.1e}")164    assert abs(e2 - SIGMA) < 2 * bar + 1e-9165    zs = (0.01, 1 / 9, 1 / 3, 3.0, 9.0, 100.0)166    val = {}167    for z in zs:168        h = [cell(m, z) for m in (64, 128, 256)]169        val[z], p = aitken(*h)170        print(f"contrast z={z:.6g} sigma={val[z]:.9f} order={p:.4f} obnosov={obnosov(z):.9f} off={val[z] - obnosov(z):+.1e}")171        assert abs(val[z] - obnosov(z)) < 6e-8172    for z in (0.01, 1 / 9, 1 / 3):173        prod = val[z] * val[1 / z if z != 1 / 3 else 3.0]174        print(f"keller z={z:.6g} sigma(z)*sigma(1/z)={prod:.9f}")175        assert abs(prod - 1) < 7e-8176    big = 1e8177    h = [cell(m, big) for m in (64, 128, 256)]178    s_inf, p = aitken(*h)179    print(f"conductor z=1e8 sigma={s_inf:.9f} order={p:.4f} sqrt3={R3:.9f} keller={e2 * s_inf:.9f}")180181def verb_side():182    for n in (11, 21, 41, 81):183        row = []184        for k in (4, 8, 16):185            row.append(n * (level(n, 1, k) - cell(k // 2, 0.0)))186        print(f"level1 N={n} N*(sigma_k(N,1)-sigma_k(inf,1)) k=4,8,16: " + " ".join(f"{x:.5f}" for x in row) + f" sigma(N,1)~{SIGMA + row[-1] / n:.6f}")187    for n in (5, 7, 9, 11, 15, 21, 31, 41):188        row = []189        for k in (1, 2, 4):190            if k * n * n > 1800 and k > 1:191                break192            row.append(level(n, 1, k) / level(n, 2, k))193        print(f"ratio N={n} sigma(N,1)/sigma(N,2) k=1,2,4: " + " ".join(f"{x:.5f}" for x in row) + f" N*(sqrt3-ratio)={n * (R3 - row[-1]):.4f}")194195def up(x):196    e = math.floor(math.log10(x))197    return f"{math.ceil(x / 10 ** e)}e{e}" if math.ceil(x / 10 ** e) < 10 else f"1e{e + 1}"198199def bounds(n):200    beta = (n + 1) / (2 * n) + (n - 1) / (n + 1)201    alpha = 2 * n / (n + 1)202    return beta, alpha203204def ratios(n, k, cap):205    g = [1.0]206    lv = 1207    while k * n ** (lv + 1) <= cap:208        if lv == 1:209            g.append(level(n, 1, k))210        g.append(level(n, lv + 1, k))211        lv += 1212    return [g[i] / g[i + 1] for i in range(1, len(g) - 1)]213214def verb_scale():215    rows = []216    for n in (3, 5, 7, 9, 11):217        beta, alpha = bounds(n)218        assert Fraction(3 * n * n + 1, 2 * n * (n + 1)) == Fraction(n + 1, 2 * n) + Fraction(n - 1, n + 1)219        seq = {k: ratios(n, k, 3200 if k == 1 else 2700) for k in (1, 2)}220        for k in (1, 2):221            print(f"scale N={n} k={k} rho(N,L) L=1..{len(seq[k])}: " + " ".join(f"{x:.5f}" for x in seq[k]))222            assert all(beta <= x <= alpha for x in seq[k])223        j = len(seq[2]) - 1224        mesh = abs(seq[1][j] - seq[2][j])225        kk = min((1, 2), key=lambda k: abs(seq[k][-1] - seq[k][-2]))226        drift = abs(seq[kk][-1] - seq[kk][-2])227        rho = seq[kk][-1]228        bar = max(mesh, drift)229        m = fill(n)230        dw = math.log(m * rho) / math.log(n)231        dbar = math.log((rho + bar) / rho) / math.log(n)232        law = 2 + math.log(3 * R3 / 4) / math.log(n)233        archie = math.log(rho) / math.log(n * n / m)234        print(f"table N={n} k={kk} rho={rho:.5f} mesh={mesh:.1e} drift={drift:.1e} bar={bar:.1e} dw_bar={dbar:.1e} bounds=[{beta:.5f}, {alpha:.5f}] dw={dw:.4f} law={law:.4f} (dw-2)logN={(dw - 2) * math.log(n):.4f} archie={archie:.4f} drift_k1={seq[1][0] - seq[1][-1]:.2e} drift_k2={seq[2][0] - seq[2][-1]:.2e} lower={(math.log(m * beta / n ** 2)):.4f} upper={(math.log(m * alpha / n ** 2)):.4f} archie_gap={math.log(3) / math.log(16 / 9) - archie:.4f} excess={(n - 1) ** 3 / (8 * n ** 3):.5f} dw-law={dw - law:.4f}")235        rows.append(f"| {n} | {kk} | " + ", ".join(f"{x:.5f}" for x in seq[kk]) + f" | `{rho:.5f} +- {up(bar)}` | `[{Fraction(3 * n * n + 1, 2 * n * (n + 1))}, {Fraction(2 * n, n + 1)}]` | `{dw:.4f} +- {up(dbar)}` | {law:.4f} |")236    print("\n".join(rows))237    print(f"limits sqrt3={R3:.5f} bounds=[1.5, 2] (dw-2)logN->{math.log(3 * R3 / 4):.4f} in [{math.log(9 / 8):.4f}, {math.log(1.5):.4f}] archie->{math.log(3) / math.log(16 / 9):.4f}")238239def verb_dim3():240    ms = (8, 16, 32, 64)241    g = [cell3(m) for m in ms]242    for m, x in zip(ms, g):243        print(f"dim3 M={m} sigma={x:.7f}")244    a1, p1 = aitken(*g[:3])245    a2, p2 = aitken(*g[1:])246    r2 = rich(g[2], g[3])247    bar = max(abs(a2 - a1), abs(a2 - r2))248    expo = math.log(1 / a2) / math.log(2)249    print(f"dim3 aitken(8,16,32)={a1:.6f} order={p1:.4f} aitken(16,32,64)={a2:.6f} order={p2:.4f} rich(32,64)={r2:.6f}")250    print(f"dim3 sigma={a2:.5f} bar={bar:.1e} exponent=log(1/sigma)/log2={expo:.4f} band=[{math.log(1 / (a2 - bar)) / math.log(2):.4f}, {math.log(1 / (a2 + bar)) / math.log(2):.4f}]")251252def verb_rate():253    ns = (9, 11, 15, 21, 31, 41)254    y = np.array([level(n, 1, 2) / level(n, 2, 2) for n in ns])255    x = np.array(ns, float)256    A = np.stack([np.ones_like(x), -1 / x, 1 / x ** (4 / 3)], axis=1)257    c, *_ = np.linalg.lstsq(A, y, rcond=None)258    res = y - A @ c259    print("rate N=" + ",".join(map(str, ns)) + " rho(N,1) k=2: " + " ".join(f"{v:.5f}" for v in y))260    print("rate N*(sqrt3-rho(N,1)): " + " ".join(f"{n * (R3 - v):.4f}" for n, v in zip(ns, y)))261    print(f"rate fit rho(N,1) = r - a/N + b/N^(4/3): r={c[0]:.4f} a={c[1]:.3f} b={c[2]:.3f} max residual={np.abs(res).max():.1e} sqrt3={R3:.4f}")262    for q in (5 / 4, 4 / 3, 3 / 2):263        A2 = np.stack([-1 / x, 1 / x ** q], axis=1)264        c2, *_ = np.linalg.lstsq(A2, y - R3, rcond=None)265        print(f"rate fit rho(N,1) = sqrt3 - a/N + b/N^{q:.4g}: a={c2[0]:.3f} b={c2[1]:.3f} max residual={np.abs(y - R3 - A2 @ c2).max():.1e}")266267def verb_dual():268    for n, lv in ((5, 1), (3, 2)):269        ext, bars = {}, {}270        for z in (0.0, 1e8):271            g = []272            for k in (8, 16, 32, 64):273                m = render(n, lv, k)274                g.append(conductance(np.where(m, 1.0, z)))275            ext[z], p = aitken(*g[1:])276            bars[z] = abs(ext[z] - aitken(*g[:3])[0])277            print(f"dual N={n} L={lv} z={z:.0e} k=8,16,32,64: " + " ".join(f"{x:.7f}" for x in g) + f" aitken={ext[z]:.7f} order={p:.3f} bar={bars[z]:.1e}")278        prod = ext[0.0] * ext[1e8]279        bar = bars[0.0] * ext[1e8] + bars[1e8] * ext[0.0]280        print(f"dual N={n} L={lv} sigma(0)*sigma(inf)={prod:.7f} bar={bar:.1e}")281        assert abs(prod - 1) < bar282283VERBS = {"cell": verb_cell, "dual": verb_dual, "side": verb_side, "rate": verb_rate, "scale": verb_scale, "dim3": verb_dim3}284285if __name__ == "__main__":286    names = sys.argv[1:] or list(VERBS)287    for name in names:288        t = time.time()289        VERBS[name]()290        print(f"{name} {time.time() - t:.1f} s", flush=True)