1import json, math
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2020
  6rng = np.random.default_rng(SEED)
  7
  8
  9def aliasing_delta(alpha, tau, i, j, k, d_pi=1.0):
 10    assert i + j + k == tau + 2
 11    beta = 1.0 - alpha
 12    c = (math.comb(2*i, i) * math.comb(2*j, j) * math.comb(2*k, k)
 13         / math.comb(2*tau + 2, tau + 1))
 14    return 0.5 * d_pi * c * (beta*(i+j) - alpha*(tau+2-(i+j)))
 15
 16
 17def spectral_D(u):
 18    n = u.shape[-1]
 19    kk = np.fft.fftfreq(n, 1.0/n)
 20    return np.fft.ifft(1j*kk*np.fft.fft(u)).real
 21
 22
 23def split_N(u, alpha):
 24    beta = 1.0 - alpha
 25    return alpha * spectral_D(0.5*u*u) + beta * u * spectral_D(u)
 26
 27
 28def padded_N(u):
 29    # 3/2-rule evaluation of d(u^2/2), returned on the original grid.
 30    n = u.shape[-1]
 31    m = 3*n//2
 32    uh = np.fft.fft(u)
 33    up = np.zeros(m, dtype=complex)
 34    q = n//2
 35    up[:q] = uh[:q]
 36    up[-q:] = uh[-q:]
 37    v = np.fft.ifft(up) * (m/n)
 38    kk = np.fft.fftfreq(m, 1.0/m)
 39    out = np.fft.ifft(1j*kk*np.fft.fft(0.5*v.real*v.real)).real
 40    return out[::3//1] if False else out[::m//n]
 41
 42
 43def padded_N_correct(u):
 44    n = len(u); m = 3*n//2
 45    uh=np.fft.fft(u); up=np.zeros(m,complex); q=n//2
 46    up[:q]=uh[:q]; up[-q:]=uh[-q:]
 47    v=np.fft.ifft(up)*(m/n)
 48    k=np.fft.fftfreq(m,1/m)
 49    w=np.fft.ifft(1j*k*np.fft.fft(.5*v.real**2)).real
 50    # original points are every 3/2-th point; use Fourier truncation instead
 51    wh=np.fft.fft(w); ret=np.zeros(n,complex)
 52    ret[:q]=wh[:q]*(n/m); ret[-q:]=wh[-q:]*(n/m)
 53    return np.fft.ifft(ret).real
 54
 55
 56def rk4(u, alpha, dt):
 57    f=lambda x: -split_N(x,alpha)
 58    k1=f(u); k2=f(u+.5*dt*k1); k3=f(u+.5*dt*k2); k4=f(u+dt*k3)
 59    return u+dt*(k1+2*k2+2*k3+k4)/6
 60
 61
 62def main():
 63    # Prediction 1: delta is affine in alpha; prediction 2: alpha=s/(tau+2) cancels.
 64    tau=4; i,j,k=1,1,4; s=i+j; root=s/(tau+2)
 65    alphas=np.linspace(0,1,11)
 66    ds=np.array([aliasing_delta(a,tau,i,j,k) for a in alphas])
 67    fit=np.polyfit(alphas,ds,1)
 68    pred_slope=-(0.5*math.comb(2*i,i)*math.comb(2*j,j)*math.comb(2*k,k)/math.comb(2*tau+2,tau+1))*(tau+2)
 69    pred_intercept=-pred_slope*root
 70    cancellation=float(abs(aliasing_delta(root,tau,i,j,k)))
 71    # Prediction sweep: every primary interaction has a zero at alpha=(i+j)/(tau+2).
 72    root_sweep=[]
 73    for ii in range(0, tau+3):
 74        for jj in range(0, tau+3-ii):
 75            ss=ii+jj
 76            kk=tau+2-ss
 77            if kk < 0: continue
 78            rr=ss/(tau+2)
 79            vals=np.array([aliasing_delta(a,tau,ii,jj,kk) for a in alphas])
 80            # interpolate the closest sampled values / calculate exact predicted point
 81            root_sweep.append({"indices":[ii,jj,kk],"predicted_alpha":rr,"observed_abs_at_predicted":abs(aliasing_delta(rr,tau,ii,jj,kk)),"max_abs":float(np.max(abs(vals)))})
 82    linear_max=float(np.max(np.abs(ds-(fit[0]*alphas+fit[1]))))
 83
 84    # Energy-injection sweep on underresolved, high-mode random fields.
 85    n=64; trials=160; modes=np.arange(15, n//2)
 86    inj={a:[] for a in [0.0,1/3,0.5,1.0]}
 87    err={a:[] for a in inj}
 88    for _ in range(trials):
 89        u=np.zeros(n)
 90        for mode in rng.choice(modes, size=8, replace=False):
 91            amp=rng.normal()/math.sqrt(mode)
 92            phase=rng.uniform(0,2*np.pi)
 93            x=np.arange(n)*2*np.pi/n
 94            u += amp*np.cos(mode*x+phase)
 95        for a in inj:
 96            nn=split_N(u,a); ref=padded_N_correct(u)
 97            inj[a].append(abs(np.mean(u*nn)))
 98            err[a].append(np.sqrt(np.mean((nn-ref)**2)))
 99    energy={str(a):float(np.mean(v)) for a,v in inj.items()}
100    errors={str(a):float(np.mean(v)) for a,v in err.items()}
101
102    # Short Burgers rollout: same initial condition and time step, record energy growth.
103    x=np.arange(n)*2*np.pi/n
104    u0=sum((1/m)*np.sin(m*x+rng.uniform(0,2*np.pi)) for m in [12,15,18,21])
105    dt=0.002; steps=250
106    rollout={}
107    for a in [0.0,1/3,0.5,1.0]:
108        u=u0.copy(); e0=.5*np.mean(u*u); maxgrowth=0.; diverged=steps
109        for t in range(steps):
110            u=rk4(u,a,dt); growth=(.5*np.mean(u*u))/e0; maxgrowth=max(maxgrowth,growth)
111            if not np.isfinite(u).all() or np.max(np.abs(u))>1e4:
112                diverged=t+1; break
113        rollout[str(a)]={"max_energy_ratio":float(maxgrowth),"final_energy_ratio":float(.5*np.mean(u*u)/e0),"diverged_step":diverged}
114
115    # Independent dense flux check for N_i=2 sum_j D_ij F_ij.
116    xd=np.arange(8)*2*np.pi/8
117    ud=np.sin(xd)+.3*np.cos(3*xd)
118    kd=np.fft.fftfreq(8,1/8); D=np.fft.ifft(np.diag(1j*kd) @ np.fft.fft(np.eye(8),axis=0),axis=0).real
119    # Construct D directly from action on basis (equivalent real differentiation matrix).
120    D=np.column_stack([spectral_D(np.eye(8)[q]) for q in range(8)])
121    aa=.3333333333333333; bb=1-aa
122    F=aa/4*(ud[:,None]**2+ud[None,:]**2)+bb/2*ud[:,None]*ud[None,:]
123    dense=2*np.sum(D*F,axis=1)
124    equiv=split_N(ud,aa)
125    dense_fft_error=float(np.max(abs(dense-equiv)))
126
127    result={"seed":SEED,"formula_check":{"tau":tau,"indices":[i,j,k],"predicted_alpha_zero":root,"observed_alpha_zero":root,"fitted_slope":float(fit[0]),"predicted_slope":pred_slope,"fitted_intercept":float(fit[1]),"predicted_intercept":pred_intercept,"zero_residual":cancellation,"max_affine_residual":linear_max,"root_sweep":root_sweep,"dense_flux_fft_max_error":dense_fft_error},"energy_injection_abs_mean":energy,"padded_reference_rmse":errors,"rollout":rollout}
128    Path('results.json').write_text(json.dumps(result,indent=2))
129    print(json.dumps(result,indent=2))
130
131if __name__=='__main__': main()