Chemical spinodal decomposition.

"""Phase-field example based on my teaching notebook PhaseField_Examples_Spectral.ipynb.
Purpose: explain spinodal decomposition and the effect of elastic strain energy.
Run: python phase_field.py [chemical|elastic] [output_directory]
Requires NumPy. CSV rows index x, columns index y; all units dimensionless.
Periodic 2D homogeneous elasticity; zero mean stress. Not a calibrated alloy.
The explicit bulk term means time-step convergence must be checked.
"""
import sys
from pathlib import Path
import numpy as np


def kernel(kx, ky, c11=4., c12=2., c44=3., eigenstrain=.2):
    """B = E:C:E - (sigma0.n).Q^-1.(sigma0.n), E=eigenstrain*I.
    Q_ik=C_ijkl*n_j*n_l; Q11=C11*nx²+C44*ny² (not C12).
    This is a 2D constitutive model, not a reduction with epsilon_zz eigenstrain.
    """
    k2 = kx*kx + ky*ky
    radius = np.sqrt(np.where(k2 == 0, 1., k2))
    nx, ny = kx/radius, ky/radius
    q11 = c11*nx*nx + c44*ny*ny
    q22 = c44*nx*nx + c11*ny*ny
    q12 = (c12+c44)*nx*ny
    det = q11*q22-q12*q12
    det = np.where(k2 == 0, 1., det)
    stress = (c11+c12)*eigenstrain
    b = 2*(c11+c12)*eigenstrain**2-stress**2*(q22*nx*nx-2*q12*nx*ny+q11*ny*ny)/det
    return np.where(k2 == 0, 0., b)


def run(mode='chemical', output='phase_output', n=None, dt=.2, steps=2000):
    out=Path(output);out.mkdir(parents=True,exist_ok=True)
    if n is None: n = 256 if mode == 'elastic' else 64
    if mode not in ('chemical','elastic'): raise ValueError('Choose chemical or elastic')
    k=2*np.pi*np.fft.fftfreq(n);kx,ky=np.meshgrid(k,k,indexing='ij');k2=kx*kx+ky*ky
    b=kernel(kx,ky) if mode == 'elastic' else np.zeros_like(k2)
    # Real-field discrete operator: symmetrize conjugate partners on even grids.
    # Only Nyquist lines can differ from the continuum directional kernel.
    partner=(-np.arange(n)) % n
    b=.5*(b+b[np.ix_(partner,partner)])
    i,j=np.meshgrid(np.arange(n),np.arange(n),indexing='ij')
    # Deterministic broadband perturbation, identical across Julia and NumPy.
    noise=np.sin((i+1)*12.9898+(j+1)*78.233)*43758.5453
    noise=.02*(noise-np.floor(noise)-.5);noise-=noise.mean();c=.5+noise
    initial_mean=c.mean();records=[]
    denominator=1+dt*k2*(k2+b)
    def record(step):
        h=np.fft.fft2(c)
        energy=np.mean(c*c*(1-c)**2)+.5*np.sum((k2+b)*np.abs(h)**2)/n**4
        records.append([step*dt,c.mean(),energy,c.min(),c.max()])
        np.savetxt(out/f'field_{step:04d}.csv',c,delimiter=',')
    record(0)
    for step in range(1,steps+1):
        h=(np.fft.fft2(c)-dt*k2*np.fft.fft2(2*c*(1-c)*(1-2*c)))/denominator
        c=np.fft.ifft2(h).real
        if step%25 == 0 or step==steps: record(step)
    assert np.isfinite(c).all()
    assert abs(c.mean()-initial_mean)<1e-10
    np.savetxt(out/'diagnostics.csv',records,delimiter=',',header='time,mean,energy_density,min,max',comments='')
    print(mode,'mean drift',abs(c.mean()-initial_mean),'energy',records[0][2],records[-1][2])
    return c,np.array(records)

if __name__=='__main__': run(*(sys.argv[1:3]))

# Phase-field example based on my teaching notebook PhaseField_Examples_Spectral.ipynb.
# Purpose: explain spinodal decomposition and the effect of elastic strain energy.
# julia phase_field.jl [chemical|elastic] [output_directory]
# Install FFTW: julia -e 'using Pkg; Pkg.add("FFTW")'
# Dimensionless periodic 2D model. CSV rows=x, columns=y.
using FFTW, Statistics, DelimitedFiles, Printf

function kernel(kx,ky; c11=4.,c12=2.,c44=3.,eigenstrain=.2)
    k2=kx .^ 2  .+  ky .^ 2
    radius=sqrt.(ifelse.(k2 .== 0,1.,k2)); nx=kx ./ radius; ny=ky ./ radius
    # Acoustic tensor Q_ik=C_ijkl*n_j*n_l: diagonal shear term is C44.
    q11=c11 .* nx .^ 2  .+  c44 .* ny .^ 2
    q22=c44 .* nx .^ 2  .+  c11 .* ny .^ 2
    q12=(c12+c44) .* nx .* ny
    det=q11 .* q22 .- q12 .^ 2; det[k2 .== 0] .= 1.
    stress=(c11+c12)*eigenstrain
    b=2*(c11+c12)*eigenstrain^2  .-  stress^2 .* (q22 .* nx .^ 2 .- 2  .* q12 .* nx .* ny .+ q11 .* ny .^ 2) ./ det
    b[k2 .== 0] .= 0.; return b
end

function run(mode="chemical",output="phase_output";n=(mode=="elastic" ? 256 : 64),dt=.2,steps=2000)
    mode in ("chemical","elastic") || error("Choose chemical or elastic")
    mkpath(output)
    k=2π .* fftfreq(n);kx=repeat(reshape(k,n,1),1,n);ky=permutedims(kx);k2=kx .^ 2 .+ ky .^ 2
    b=mode=="elastic" ? kernel(kx,ky) : zeros(n,n)
    # Conjugate-pair symmetrization of the discrete real-field kernel.
    partner=mod.(-(0:n-1),n) .+ 1
    b=0.5 .* (b .+ b[partner,partner])
    noise=[sin(i*12.9898+j*78.233)*43758.5453 for i=1:n,j=1:n]
    noise=.02 .* (noise .- floor.(noise) .-  .5);noise .-=mean(noise);c=.5 .+ noise
    initial_mean=mean(c);records=Vector{Float64}[]
    denominator=1  .+ dt .* k2 .* (k2 .+ b)
    function record(step)
        h=fft(c);energy=mean(c .^ 2 .* (1  .- c) .^ 2)+.5*sum((k2 .+ b) .* abs2.(h))/n^4
        push!(records,[step*dt,mean(c),energy,minimum(c),maximum(c)])
        writedlm(joinpath(output,@sprintf("field_%04d.csv",step)),c,',')
    end
    record(0)
    for step=1:steps
        h=(fft(c) .- dt .* k2 .* fft(2  .* c .* (1  .- c) .* (1  .- 2  .* c))) ./ denominator
        c=real.(ifft(h))
        (step%25==0 || step==steps) && record(step)
    end
    @assert all(isfinite,c)
    @assert abs(mean(c)-initial_mean)<1e-10
    open(joinpath(output,"diagnostics.csv"),"w") do io
        println(io,"time,mean,energy_density,min,max");writedlm(io,permutedims(hcat(records...)),',')
    end
    println(mode," mean drift ",abs(mean(c)-initial_mean)," energy ",records[1][3]," → ",records[end][3])
    return c,records
end
if abspath(PROGRAM_FILE) == @__FILE__
    run(ARGS...)
end