The file src/32-dark/flux_exact.py of the validation package of the paper 32-dark, rendered with an anchor per line; the file itself is beside this page.
1# Construction 1, EXACT no-go (framed-rational, zero floating point).2# Linear Kuramoto statics L u = source on a 5^3 box, Dirichlet boundary, exact rationals.3# Verifies the discrete Gauss law: net synchronisation flux through any closed surface4# enclosing the source = enclosed cardinality, EXACTLY; a source-free (noise) field carries5# zero net flux; hence source+noise has the same enclosed flux as the source. The mean force6# is therefore independent of the floor noise: Lemma fluxnoise, exact.7from fractions import Fraction as Fr8import itertools, random9N=510pts=list(itertools.product(range(N),repeat=3))11def interior(p): return all(0<c<N-1 for c in p)12I=[p for p in pts if interior(p)]; pos={p:i for i,p in enumerate(I)}; M=len(I)13ctr=(2,2,2)14def nbrs(p):15 x,y,z=p16 for d in((1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)):17 yield (x+d[0],y+d[1],z+d[2])18def solve(A,b):19 n=len(b); A=[r[:] for r in A]; b=b[:]20 for c in range(n):21 piv=next(r for r in range(c,n) if A[r][c]!=0)22 A[c],A[piv]=A[piv],A[c]; b[c],b[piv]=b[piv],b[c]23 inv=Fr(1)/A[c][c]; A[c]=[v*inv for v in A[c]]; b[c]*=inv24 for r in range(n):25 if r!=c and A[r][c]!=0:26 f=A[r][c]; A[r]=[a-f*x for a,x in zip(A[r],A[c])]; b[r]-=f*b[c]27 return b28def field(rho, bvals): # bvals: dict of fixed boundary values29 A=[[Fr(0)]*M for _ in range(M)]; b=[Fr(0)]*M30 for p in I:31 i=pos[p]; A[i][i]=Fr(6)32 for q in nbrs(p):33 if q in pos: A[i][pos[q]]-=134 else: b[i]+=bvals.get(q,Fr(0))35 b[i]+=rho.get(p,Fr(0))36 x=solve(A,b); U=dict(bvals)37 for p in I: U[p]=x[pos[p]]38 return U39def flux(U,R):40 tot=Fr(0)41 for p in pts:42 if max(abs(p[0]-2),abs(p[1]-2),abs(p[2]-2))>R: continue43 for q in nbrs(p):44 if max(abs(q[0]-2),abs(q[1]-2),abs(q[2]-2))>R:45 tot += U.get(p,Fr(0))-U.get(q,Fr(0))46 return tot47zero_b={p:Fr(0) for p in pts if not interior(p)}48Us=field({ctr:Fr(1)}, zero_b)49print("EXACT discrete Gauss law (rational):")50for R in (0,1): print(f" unit point source: flux(R={R}) = {flux(Us,R)} (enclosed = 1)")51random.seed(1)52nb_b={p:Fr(random.randint(-9,9),4) for p in pts if not interior(p)}53Un=field({}, nb_b)54print(f" source-free noise field: flux(R=1) = {flux(Un,1)} (no enclosed source -> 0)")55Usn=field({ctr:Fr(1)}, nb_b)56print(f" source+noise: flux(R=1) = {flux(Usn,1)} == source-alone {flux(Us,1)} ? {flux(Usn,1)==flux(Us,1)}")57print("PASS" if flux(Us,1)==Fr(1) and flux(Un,1)==Fr(0) and flux(Usn,1)==Fr(1) else "FAIL")58