1import json, math, time
 2import numpy as np
 3import torch
 4from hodge_dual import curl_2d, divergence_2d
 5
 6torch.set_default_dtype(torch.float64)
 7torch.manual_seed(7); np.random.seed(7)
 8
 9def spectrum(N, eps):
10    h=1.0/N
11    # central-difference curl symbol; zero mode is gauge nullspace
12    vals=[]
13    for k in range(N):
14        sx=math.sin(2*math.pi*k/N)/h
15        for l in range(N):
16            sy=math.sin(2*math.pi*l/N)/h
17            vals.append((sx*sx+sy*sy)/eps)
18    return np.array(vals)/(N*N)
19
20def divergence_check():
21    N=24; h=1/N
22    A=torch.randn(1,N,N,3)*3.7
23    p0=torch.randn(1,N,N,2)
24    r=divergence_2d(curl_2d(A,h),h)
25    # random amplitudes demonstrate independence from potential magnitude
26    rows=[]
27    for scale in [0.,1e-4,1.,1e4]:
28        rr=divergence_2d(curl_2d(A*scale,h),h)
29        rows.append(float(rr.abs().max()))
30    return rows
31
32def gd_mode(N, eps, eta, steps=80):
33    h=1/N
34    # one Fourier-like real perturbation, with p0=0; energy is purely quadratic
35    x=torch.arange(N).reshape(1,N,1).double(); y=torch.arange(N).reshape(1,1,N).double()
36    A=torch.zeros(1,N,N,3)
37    A[...,2]=torch.sin(2*math.pi*6*x/N)*torch.sin(2*math.pi*6*y/N)
38    initial=0.5*(curl_2d(A,h)**2).sum()/eps/(N*N)
39    vals=[]
40    for _ in range(steps):
41        A.requires_grad_(True)
42        loss=0.5*(curl_2d(A,h)**2).sum()/eps/(N*N)
43        g=torch.autograd.grad(loss,A)[0]
44        with torch.no_grad(): A=A-eta*g
45        vals.append(float(loss))
46    return float(initial), vals[-1], vals
47
48def primal_compare(N=24):
49    # Same central differences: primal Poisson objective and dual quadratic have
50    # identical positive conditioning scale (dual additionally enforces div p).
51    vals=spectrum(N,1.0); positive=vals[vals>1e-12]
52    return {"primal_condition_number":float(positive.max()/positive.min()),
53            "dual_condition_number":float(positive.max()/positive.min()),
54            "dual_gauss_residual":max(divergence_check())}
55
56def main():
57    N=24
58    divs=divergence_check()
59    # Prediction 1: div(curl A)=0 to roundoff at every amplitude.
60    div_pred=0.0
61    # Prediction 2: lambda_max(eps)*eps is constant.
62    epss=[0.5,1.,2.,4.]
63    lams=[float(spectrum(N,e).max()) for e in epss]
64    products=[e*l for e,l in zip(epss,lams)]
65    # Prediction 3: GD changes from contraction to divergence at eta*lambda_max=2.
66    lam=lams[1]
67    eta_crit=2/lam
68    tests=[]
69    for mult in [0.90,0.99,1.01,1.10]:
70        initial, final, _=gd_mode(N,1.,mult*eta_crit,steps=30)
71        tests.append({"eta_over_critical":mult,"final_over_initial":final/initial,
72                      "stable_observed":bool(final < initial*10)})
73    out={"grid":N,"divergence_max_abs_by_amplitude":divs,
74         "prediction_1":{"predicted":"0 exactly (up to floating point)","observed_max":max(divs)},
75         "prediction_2":{"predicted":"lambda_max * eps is constant","eps":epss,
76                          "lambda_max":lams,"products":products,
77                          "relative_spread":(max(products)-min(products))/np.mean(products)},
78         "prediction_3":{"predicted":"transition at eta*lambda_max=2","lambda_max":lam,
79                          "eta_critical":eta_crit,"tests":tests},
80         "baseline_context":primal_compare(N),
81         "note":"The primal and dual constant-coefficient quadratic operators have the same nonzero Fourier condition number; the dual advantage tested here is exact constraint satisfaction, not a universal conditioning improvement."}
82    with open('results.json','w') as f: json.dump(out,f,indent=2)
83    print(json.dumps(out,indent=2))
84
85if __name__=='__main__': main()