rdeasy.py

import math
import numpy as np

N = 16
gridsize = N*N*N
t = 0
t0 = 0
tf = 10
dt = 0.05
L = 20 # box size
dx = L / N
lambad_ = 9*math.pow(10,-14) # self coupling inflton
gl = 200 #  resonance parameter g^2/lamba_


phi = np.ones((N, N, N))/1000
dphi = np.zeros((N, N, N))
chi = np.zeros((N, N, N))
dchi = np.zeros((N,N,N))

def inc(k):
    if k == N - 1:
        return 0
    return k + 1


def dec(k):
    if k == N - 1:
        return 0
    return k+1


def laplace_b(f, x, y, z):
    return f[x][y][inc(z)] + f[x][y][dec(z)]+ f[x][inc(y)][z] + f[x][dec(y)][z] +f[inc(x)][y][inc(z)] + f[dec(x)][y][z]    - 6 * f[x][y][z]


def laplace(f, x, y, z):
    # print("laplace"+str(x)+","+str(y)+","+str(z))
    return f[x][y][z+1] + f[x][y][z] -f[x][y+1][z] + f[x][y-1][z] + f[x+1][y][z] + f[x-1][y][z] - 6 * f[x][y][z]


def dVdphi(phi_v, dphi_v, chi_v, dchi_v):
    # print("phi_v+",phi_v)
    return math.pow(phi_v, 2) + math.pow(chi_v, 2) * phi_v

def dVdchi(phi_v, dphi_v, chi_v, dchi_v):
    # print("chi_v+", phi_v)
    return gl*math.pow(phi_v,2)*chi_v   


def evolve_fields(f,df,dt):
    f += f * df
    return f

def evolve_d_phi(phi,dphi,chi,dchi,x,y,z,dt):
    return dt * laplace(phi, x, y, z) + dVdphi(phi[x][y][z], dphi[x][y][z], chi[x][y][z], dchi[x][y][z])
     
def evolve_d_chi(phi, dphi, chi, dchi, x, y, z, dt):
    return dt * laplace(phi, x, y, z) + dVdchi(phi[x][y][z], dphi[x][y][z], chi[x][y][z], dchi[x][y][z])


def evolve_d_filed(phi, dphi, chi, dchi, dt):
    for x in range(1,N-1):
        for y in range(1,N-1):
            for z in range(1,N-1):
                dphi[x, y, z] = evolve_d_phi(phi, dphi, chi, dchi, x, y, z,dt)
                dchi[x, y, z] = evolve_d_chi(
                    phi, dphi, chi, dchi, x, y, z, dt)

# 端っこの閭里todo
 

def evolve_main(t0, tf, phi, dphi, chi, dchi, dt, gridsize):
    t = t0
    evolve_d_filed(phi, dchi, chi, dchi, 0.5*dt)
    while (t <= tf):
        evolve_d_filed(phi, dphi, chi, dchi, dt)
        phi = evolve_fields(phi, dphi,dt)
        chi = evolve_fields(chi, dchi, dt)
        t += dt
        recordData(t0, tf, phi, dphi, chi, dchi, dt, gridsize)


def fieldTotal(f):
    r = 0
    for x in range(N):
        for y in range(N):
            for z in range(N):
                r = f[x][y][z]
    
    return r


def fieldTotalPow2(f,minus):
    r = 0
    for x in range(N):
        for y in range(N):
            for z in range(N):
                r = f[x][y][z]*f[x][y][z] - minus

    return r


def fieldTotalProd(f, df, minus):
    r = 0
    for x in range(N):
        for y in range(N):
            for z in range(N):
                r = f[x][y][z]*df[x][y][z] - minus

    return r

def meansvars(f, gridsize):
    av = fieldTotal(f) / gridsize
    av_sq = av * av
    
    var = fieldTotalPow2(f, av_sq)
    return av, var

def gradientEnergy(f, df, gridsize, dx):
    # print(f)
    gradient = 0
    for x in range(1,N-1):
        for y in range(1,N-1):
            for z in range(1, N - 1):
                # print(""+str(x)+","+str(y)+","+str(z))
                gradient -= f[x][y][z] * laplace(f, x, y, z)

    for i in range(N):
        for j in range(N):
            gradient -= f[i][j][0] * laplace_b(f, i, j, 0) + \
                f[i][j][N-1]*laplace_b(f, i, j, N-1)
            if j == 0 or j == N - 1:
                continue
            gradient -= f[i][0][j] * laplace_b(
                f, i, 0, j) + f[i][N-1][j]*laplace_b(f, i, N-1, j)
            if i == 0 or i == N - 1:
                continue
            gradient -= f[0][i][j] * laplace_b(
                f, 0, i, j) + f[N-1][i][j]*laplace_b(f, N-1, i, j)

    norm = 1/math.pow(dx,2)
    return 0.5 * gradient*norm/gridsize
                




def energy(f, df, gridsize,dx):
    var = fieldTotalPow2(f, 0) / gridsize
    vard = fieldTotalPow2(df, 0) / gridsize
    ffd = fieldTotalProd(f,df,0) / gridsize
    deriv_energy = 0.5 * vard  # model
    grad_energy = gradientEnergy(f,df,gridsize,dx)


phi_meanes = []
phi_var = []
phi_energy = []
chi_means = []
chi_var = []
chi_energy = []


def recordData(t0, tf, phi, dphi, chi, dchi, dt, gridsize):
    if t>0:
        evolve_fields(-0.5 * dt)
    
    means_, vars_ = meansvars(phi, gridsize)
    energy_ = energy(phi, dphi, gridsize, dx)
    phi_meanes.append(means_)
    phi_var.append(vars_)
    phi_energy.append(energy_)

    means_, vars_ = meansvars(chi, gridsize)
    energy_ = energy(chi, dchi,gridsize, dx)
    chi_means.append(means_)
    chi_var.append(vars_)
    chi_energy.append(energy_)


    if t > 0:
        evolve_fields(0.5 * dt)
        
    return means_, vars_, energy_


def initialize():
    recordData(t0, tf, phi, dphi, chi, dchi, dt, gridsize)
    return


def main():
    initialize()
    evolve_main(t0, tf, phi, dphi, chi, dchi, dt, gridsize)

    print("phi_means")
    print(phi_meanes)

    print("phi_vars")
    print(phi_vars)

    print("phi_energy")
    print(phi_energy)

    print("chi_means")
    print(chi_means)

    print("chi_vars")
    print(chi_vars)

    print("chi_energy")
    print(chi_energy)

main()