lattice/sfcf.jl
1
0
Fork 0
mirror of https://igit.ific.uv.es/fernando.p.csic.es/sfcf.jl.git synced 2025-11-30 01:46:52 +01:00
Code Issues Projects Releases Packages Wiki Activity

Split of the code in src

Browse source
This commit is contained in:
Fernando P. Panadero 2024-02-21 17:05:34 +01:00
commit 3cb578d033
4 changed files with 459 additions and 349 deletions
Download patch file Download diff file Expand all files Collapse all files

358
sfcf.jl
View file

@ -4,360 +4,20 @@ using TimerOutputs
using ArgParse
using CUDA
function parse_commandline()
s = ArgParseSettings()
include("./src/io.jl")
include("./src/props.jl")
include("./src/corr.jl")
@add_arg_table s begin
"-i"
help = "Input parameters file"
required = true
arg_type = String
"-c"
help = "Gauge configuration file"
required = true
arg_type = String
"--cern"
help = "Config written with the export_cnfg_cern() convention"
action = :store_true
end
return parse_args(s)
end
parsed_args = parse_commandline()
println("--------------------------------------------------")
println("Reading input file from:", parsed_args["i"], "...")
params = TOML.parsefile(parsed_args["i"])
lp = SpaceParm{4}(tuple(params["Space"]["size"]...), tuple(params["Space"]["blocks"]...),BC_SF_ORBI, (0,0,0,0,0,0))
gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(params["Space"]["cG"],0.0),params["Space"]["phiT"],lp.iL);
dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,params["Fermion"]["csw"],ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,params["Fermion"]["ct"]);
dws = DiracWorkspace(SU3fund,Float64,lp);
ymws = YMworkspace(SU3,Float64,lp);
println("Reading gauge field from: ", parsed_args["c"], "...")
if !parsed_args["cern"]
U,_ = read_cnfg(parsed_args["c"])
else
U = vector_field(SU3{Float64}, lp);
file = open(parsed_args["c"])
for t in 1:lp.iL[4]
for i in 1:lp.iL[1]
for j in 1:lp.iL[2]
for k in 1:lp.iL[3]
for d in [4,1,2,3]
f,r = point_index(CartesianIndex((i,j,k,t)),lp)
#a11 !!
re11 = read(file,Float64)
co11 = read(file,Float64)
#a12 !!
re12 = read(file,Float64)
co12 = read(file,Float64)
#a13 !!
re13 = read(file,Float64)
co13 = read(file,Float64)
#a21 !!
re21 = read(file,Float64)
co21 = read(file,Float64)
#a22 !!
re22 = read(file,Float64)
co22 = read(file,Float64)
#a23 !!
re23 = read(file,Float64)
co23 = read(file,Float64)
#a31
re31 = read(file,Float64)
co31 = read(file,Float64)
#a32
re32 = read(file,Float64)
co32 = read(file,Float64)
#a33
re33 = read(file,Float64)
co33 = read(file,Float64)
CUDA.@allowscalar (U[f,d,r] = SU3{Float64}(re11 + im*co11, re12 + im*co12, re13 + im*co13,
re21 + im*co21, re22 + im*co22, re23 + im*co23))
end
end
end
end
end
length(read(file)) == (prod(lp.iL[1:3])*4*8*9*2) ? nothing : error("File not fully read")
close(file)
@timeit "Input reading" begin
read_input()
load_structs()
U = load_gauge_field()
end
Csw!(dws, U, gp, lp)
println("--------------------------------------------------")
@timeit "Propagator computation" compute_propagators()
println("Parameters:")
println("Lattice size: ", lp.iL)
println("Phi0 = ", params["Space"]["phi0"])
println("PhiT = ", params["Space"]["phiT"])
println("cG = ", gp.cG[1])
println("kappa = ", params["Fermion"]["kappa"])
println("theta = ", params["Fermion"]["theta"])
println("csw = ", dpar.csw)
println("ct = ", dpar.ct)
println("tolerance = ", params["Solver"]["tolerance"])
println("maxiter = ", params["Solver"]["maxiter"])
println("--------------------------------------------------")
psi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
println("Computing propagators from t=0 to the bulk")
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A11 = Array(psi)
C11 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
println("CG converged for c=2 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A21 = Array(psi)
C21 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 1)
println("CG converged for c=3 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A31 = Array(psi)
C31 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 2)
println("CG converged for c=1 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A12 = Array(psi)
C12 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 2)
println("CG converged for c=2 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A22 = Array(psi)
C22 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 2)
println("CG converged for c=3 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
A32 = Array(psi)
C32 = bndtobnd(psi, U, dpar, dws, lp)
println("Computing propagators from t=T to the bulk")
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B11 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
println("CG converged for c=2 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B21 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 1)
println("CG converged for c=3 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B31 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 2)
println("CG converged for c=1 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B12 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 2)
println("CG converged for c=2 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B22 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 2)
println("CG converged for c=3 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
B32 = Array(psi)
println("--------------------------------------------------")
println("Computing correlators")
println("Computing f1...")
file = open("./output/"*basename(parsed_args["c"])*".f1","w+")
f1 = (norm2(C11) + norm2(C21) + norm2(C31) + norm2(C12) + norm2(C22) + norm2(C32))
println(file," ",f1)
flush(file)
close(file)
println("Computing k1...")
file = open("./output/"*basename(parsed_args["c"])*".k1","w+")
k1 = (dot(C11,dmul(Gamma{6},imm(C12))) + dot(C11,dmul(Gamma{7},-C12)) + dot(C11,dmul(Gamma{8}, imm(C11)))
+dot(C12,dmul(Gamma{6},imm(C11))) + dot(C12,dmul(Gamma{7}, C11)) + dot(C12,dmul(Gamma{8},mimm(C12)))
+dot(C21,dmul(Gamma{6},imm(C22))) + dot(C21,dmul(Gamma{7},-C22)) + dot(C21,dmul(Gamma{8}, imm(C21)))
+dot(C22,dmul(Gamma{6},imm(C21))) + dot(C22,dmul(Gamma{7}, C21)) + dot(C22,dmul(Gamma{8},mimm(C22)))
+dot(C31,dmul(Gamma{6},imm(C32))) + dot(C31,dmul(Gamma{7},-C32)) + dot(C31,dmul(Gamma{8}, imm(C31)))
+dot(C32,dmul(Gamma{6},imm(C31))) + dot(C32,dmul(Gamma{7}, C31)) + dot(C32,dmul(Gamma{8},mimm(C32))) )/3
println(file," ",k1)
flush(file)
close(file)
println("Computing fP...")
file = open("./output/"*basename(parsed_args["c"])*".fP","w+")
fP = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
fP[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
fP[t] += (norm2(A11[b,r]) + norm2(A21[b,r]) + norm2(A31[b,r]) + norm2(A12[b,r]) + norm2(A22[b,r]) + norm2(A32[b,r]))/prod(lp.iL[1:3])
end end end
println(file,t," ",fP[t])
end
flush(file)
close(file)
println("Computing fA...")
file = open("./output/"*basename(parsed_args["c"])*".fA","w+")
fA = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
fA[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
fA[t] -= (dot(A11[b,r],dmul(Gamma{4},A11[b,r])) + dot(A21[b,r],dmul(Gamma{4},A21[b,r])) + dot(A31[b,r],dmul(Gamma{4},A31[b,r]))
+ dot(A12[b,r],dmul(Gamma{4},A12[b,r])) + dot(A22[b,r],dmul(Gamma{4},A22[b,r])) + dot(A32[b,r],dmul(Gamma{4},A32[b,r])))/prod(lp.iL[1:3])
end end end
println(file,t," ",fA[t])
end
flush(file)
close(file)
println("Computing gP...")
file = open("./output/"*basename(parsed_args["c"])*".gP","w+")
gP = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
gP[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
gP[t] += (norm2(B11[b,r]) + norm2(B21[b,r]) + norm2(B31[b,r]) + norm2(B12[b,r]) + norm2(B22[b,r]) + norm2(B32[b,r]))/prod(lp.iL[1:3])
end end end
println(file,t," ",gP[t])
end
flush(file)
close(file)
println("Computing gA...")
file = open("./output/"*basename(parsed_args["c"])*".gA","w+")
gA = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
gA[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
gA[t] += (dot(B11[b,r],dmul(Gamma{4},B11[b,r])) + dot(B21[b,r],dmul(Gamma{4},B21[b,r])) + dot(B31[b,r],dmul(Gamma{4},B31[b,r]))
+ dot(B12[b,r],dmul(Gamma{4},B12[b,r])) + dot(B22[b,r],dmul(Gamma{4},B22[b,r])) + dot(B32[b,r],dmul(Gamma{4},B32[b,r])))/prod(lp.iL[1:3])
end end end
println(file,t," ",gA[t])
end
flush(file)
close(file)
#### print f1
println("Computing kV...")
file = open("./output/"*basename(parsed_args["c"])*".kV","w+")
kv = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
kv[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
kv[t] += ( dot(A11[b,r],dmul(Gamma{6},imm(A12[b,r]))) + dot(A11[b,r],dmul(Gamma{7},-A12[b,r])) + dot(A11[b,r],dmul(Gamma{8}, imm(A11[b,r])))
+dot(A12[b,r],dmul(Gamma{6},imm(A11[b,r]))) + dot(A12[b,r],dmul(Gamma{7}, A11[b,r])) + dot(A12[b,r],dmul(Gamma{8},mimm(A12[b,r])))
+dot(A21[b,r],dmul(Gamma{6},imm(A22[b,r]))) + dot(A21[b,r],dmul(Gamma{7},-A22[b,r])) + dot(A21[b,r],dmul(Gamma{8}, imm(A21[b,r])))
+dot(A22[b,r],dmul(Gamma{6},imm(A21[b,r]))) + dot(A22[b,r],dmul(Gamma{7}, A21[b,r])) + dot(A22[b,r],dmul(Gamma{8},mimm(A22[b,r])))
+dot(A31[b,r],dmul(Gamma{6},imm(A32[b,r]))) + dot(A31[b,r],dmul(Gamma{7},-A32[b,r])) + dot(A31[b,r],dmul(Gamma{8}, imm(A31[b,r])))
+dot(A32[b,r],dmul(Gamma{6},imm(A31[b,r]))) + dot(A32[b,r],dmul(Gamma{7}, A31[b,r])) + dot(A32[b,r],dmul(Gamma{8},mimm(A32[b,r]))) )/(3*prod(lp.iL[1:3]))
end end end
println(file,t," ",kv[t])
end
flush(file)
close(file)
println("Computing lV...")
file = open("./output/"*basename(parsed_args["c"])*".lV","w+")
lv = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
lv[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
lv[t] += ( dot(B11[b,r],dmul(Gamma{6},mimm(B12[b,r]))) + dot(B11[b,r],dmul(Gamma{7}, B12[b,r])) + dot(B11[b,r],dmul(Gamma{8},mimm(B11[b,r])))
+dot(B12[b,r],dmul(Gamma{6},mimm(B11[b,r]))) + dot(B12[b,r],dmul(Gamma{7},-B11[b,r])) + dot(B12[b,r],dmul(Gamma{8}, imm(B12[b,r])))
+dot(B21[b,r],dmul(Gamma{6},mimm(B22[b,r]))) + dot(B21[b,r],dmul(Gamma{7}, B22[b,r])) + dot(B21[b,r],dmul(Gamma{8},mimm(B21[b,r])))
+dot(B22[b,r],dmul(Gamma{6},mimm(B21[b,r]))) + dot(B22[b,r],dmul(Gamma{7},-B21[b,r])) + dot(B22[b,r],dmul(Gamma{8}, imm(B22[b,r])))
+dot(B31[b,r],dmul(Gamma{6},mimm(B32[b,r]))) + dot(B31[b,r],dmul(Gamma{7}, B32[b,r])) + dot(B31[b,r],dmul(Gamma{8},mimm(B31[b,r])))
+dot(B32[b,r],dmul(Gamma{6},mimm(B31[b,r]))) + dot(B32[b,r],dmul(Gamma{7},-B31[b,r])) + dot(B32[b,r],dmul(Gamma{8}, imm(B32[b,r]))) )/(3*prod(lp.iL[1:3]))
end end end
println(file,t," ",lv[t])
end
flush(file)
close(file)
println("Computing kT...")
file = open("./output/"*basename(parsed_args["c"])*".kT","w+")
kt = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
kt[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
kt[t] +=im*(dot(A11[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A12[b,r])))) + dot(A11[b,r],dmul(Gamma{11},dmul(Gamma{5},-A12[b,r]))) + dot(A11[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A11[b,r]))))
+dot(A12[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A11[b,r])))) + dot(A12[b,r],dmul(Gamma{11},dmul(Gamma{5}, A11[b,r]))) + dot(A12[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A12[b,r]))))
+dot(A21[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A22[b,r])))) + dot(A21[b,r],dmul(Gamma{11},dmul(Gamma{5},-A22[b,r]))) + dot(A21[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A21[b,r]))))
+dot(A22[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A21[b,r])))) + dot(A22[b,r],dmul(Gamma{11},dmul(Gamma{5}, A21[b,r]))) + dot(A22[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A22[b,r]))))
+dot(A31[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A32[b,r])))) + dot(A31[b,r],dmul(Gamma{11},dmul(Gamma{5},-A32[b,r]))) + dot(A31[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A31[b,r]))))
+dot(A32[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A31[b,r])))) + dot(A32[b,r],dmul(Gamma{11},dmul(Gamma{5}, A31[b,r]))) + dot(A32[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A32[b,r])))) )/(3*prod(lp.iL[1:3]))
end end end
println(file,t," ",kt[t])
end
flush(file)
close(file)
println("Computing lT...")
file = open("./output/"*basename(parsed_args["c"])*".lT","w+")
lt = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
lt[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
lt[t] +=-im*(dot(B11[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B12[b,r])))) + dot(B11[b,r],dmul(Gamma{11},dmul(Gamma{5}, B12[b,r]))) + dot(B11[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B11[b,r]))))
+dot(B12[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B11[b,r])))) + dot(B12[b,r],dmul(Gamma{11},dmul(Gamma{5},-B11[b,r]))) + dot(B12[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B12[b,r]))))
+dot(B21[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B22[b,r])))) + dot(B21[b,r],dmul(Gamma{11},dmul(Gamma{5}, B22[b,r]))) + dot(B21[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B21[b,r]))))
+dot(B22[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B21[b,r])))) + dot(B22[b,r],dmul(Gamma{11},dmul(Gamma{5},-B21[b,r]))) + dot(B22[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B22[b,r]))))
+dot(B31[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B32[b,r])))) + dot(B31[b,r],dmul(Gamma{11},dmul(Gamma{5}, B32[b,r]))) + dot(B31[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B31[b,r]))))
+dot(B32[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B31[b,r])))) + dot(B32[b,r],dmul(Gamma{11},dmul(Gamma{5},-B31[b,r]))) + dot(B32[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B32[b,r])))) )/(3*prod(lp.iL[1:3]))
end end end
println(file,t," ",lt[t])
end
flush(file)
close(file)
@timeit "Correlators contractions" compute_correlators()
print("\n\n")

213
src/corr.jl Normal file
View file

@ -0,0 +1,213 @@
using LatticeGPU
using TOML
using TimerOutputs
using ArgParse
using CUDA
using BenchmarkTools
function fP_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
fP[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] += (norm2(A11[b,r]) + norm2(A21[b,r]) + norm2(A31[b,r]) + norm2(A12[b,r]) + norm2(A22[b,r]) + norm2(A32[b,r]))/prod(lp.iL[1:3])
end end end
end
return cf
end
function fA_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] -= (dot(A11[b,r],dmul(Gamma{4},A11[b,r])) + dot(A21[b,r],dmul(Gamma{4},A21[b,r])) + dot(A31[b,r],dmul(Gamma{4},A31[b,r]))
+ dot(A12[b,r],dmul(Gamma{4},A12[b,r])) + dot(A22[b,r],dmul(Gamma{4},A22[b,r])) + dot(A32[b,r],dmul(Gamma{4},A32[b,r])))/prod(lp.iL[1:3])
end end end
end
return cf
end
function gP_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] += (norm2(B11[b,r]) + norm2(B21[b,r]) + norm2(B31[b,r]) + norm2(B12[b,r]) + norm2(B22[b,r]) + norm2(B32[b,r]))/prod(lp.iL[1:3])
end end end
end
return cf
end
function ga_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] += (dot(B11[b,r],dmul(Gamma{4},B11[b,r])) + dot(B21[b,r],dmul(Gamma{4},B21[b,r])) + dot(B31[b,r],dmul(Gamma{4},B31[b,r]))
+ dot(B12[b,r],dmul(Gamma{4},B12[b,r])) + dot(B22[b,r],dmul(Gamma{4},B22[b,r])) + dot(B32[b,r],dmul(Gamma{4},B32[b,r])))/prod(lp.iL[1:3])
end end end
end
return cf
end
function kv_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] += ( dot(A11[b,r],dmul(Gamma{6},imm(A12[b,r]))) + dot(A11[b,r],dmul(Gamma{7},-A12[b,r])) + dot(A11[b,r],dmul(Gamma{8}, imm(A11[b,r])))
+dot(A12[b,r],dmul(Gamma{6},imm(A11[b,r]))) + dot(A12[b,r],dmul(Gamma{7}, A11[b,r])) + dot(A12[b,r],dmul(Gamma{8},mimm(A12[b,r])))
+dot(A21[b,r],dmul(Gamma{6},imm(A22[b,r]))) + dot(A21[b,r],dmul(Gamma{7},-A22[b,r])) + dot(A21[b,r],dmul(Gamma{8}, imm(A21[b,r])))
+dot(A22[b,r],dmul(Gamma{6},imm(A21[b,r]))) + dot(A22[b,r],dmul(Gamma{7}, A21[b,r])) + dot(A22[b,r],dmul(Gamma{8},mimm(A22[b,r])))
+dot(A31[b,r],dmul(Gamma{6},imm(A32[b,r]))) + dot(A31[b,r],dmul(Gamma{7},-A32[b,r])) + dot(A31[b,r],dmul(Gamma{8}, imm(A31[b,r])))
+dot(A32[b,r],dmul(Gamma{6},imm(A31[b,r]))) + dot(A32[b,r],dmul(Gamma{7}, A31[b,r])) + dot(A32[b,r],dmul(Gamma{8},mimm(A32[b,r]))) )/(3*prod(lp.iL[1:3]))
end end end
end
return cf
end
function lv_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] += ( dot(B11[b,r],dmul(Gamma{6},mimm(B12[b,r]))) + dot(B11[b,r],dmul(Gamma{7}, B12[b,r])) + dot(B11[b,r],dmul(Gamma{8},mimm(B11[b,r])))
+dot(B12[b,r],dmul(Gamma{6},mimm(B11[b,r]))) + dot(B12[b,r],dmul(Gamma{7},-B11[b,r])) + dot(B12[b,r],dmul(Gamma{8}, imm(B12[b,r])))
+dot(B21[b,r],dmul(Gamma{6},mimm(B22[b,r]))) + dot(B21[b,r],dmul(Gamma{7}, B22[b,r])) + dot(B21[b,r],dmul(Gamma{8},mimm(B21[b,r])))
+dot(B22[b,r],dmul(Gamma{6},mimm(B21[b,r]))) + dot(B22[b,r],dmul(Gamma{7},-B21[b,r])) + dot(B22[b,r],dmul(Gamma{8}, imm(B22[b,r])))
+dot(B31[b,r],dmul(Gamma{6},mimm(B32[b,r]))) + dot(B31[b,r],dmul(Gamma{7}, B32[b,r])) + dot(B31[b,r],dmul(Gamma{8},mimm(B31[b,r])))
+dot(B32[b,r],dmul(Gamma{6},mimm(B31[b,r]))) + dot(B32[b,r],dmul(Gamma{7},-B31[b,r])) + dot(B32[b,r],dmul(Gamma{8}, imm(B32[b,r]))) )/(3*prod(lp.iL[1:3]))
end end end
end
return cf
end
function kt_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] +=im*(dot(A11[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A12[b,r])))) + dot(A11[b,r],dmul(Gamma{11},dmul(Gamma{5},-A12[b,r]))) + dot(A11[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A11[b,r]))))
+dot(A12[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A11[b,r])))) + dot(A12[b,r],dmul(Gamma{11},dmul(Gamma{5}, A11[b,r]))) + dot(A12[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A12[b,r]))))
+dot(A21[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A22[b,r])))) + dot(A21[b,r],dmul(Gamma{11},dmul(Gamma{5},-A22[b,r]))) + dot(A21[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A21[b,r]))))
+dot(A22[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A21[b,r])))) + dot(A22[b,r],dmul(Gamma{11},dmul(Gamma{5}, A21[b,r]))) + dot(A22[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A22[b,r]))))
+dot(A31[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A32[b,r])))) + dot(A31[b,r],dmul(Gamma{11},dmul(Gamma{5},-A32[b,r]))) + dot(A31[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A31[b,r]))))
+dot(A32[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A31[b,r])))) + dot(A32[b,r],dmul(Gamma{11},dmul(Gamma{5}, A31[b,r]))) + dot(A32[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A32[b,r])))) )/(3*prod(lp.iL[1:3]))
end end end
end
return cf
end
function lt_fun()
cf = Vector{Complex{Float64}}(undef,lp.iL[4])
for t in 1:lp.iL[4]
cf[t] = 0.0
for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
cf[t] +=-im*(dot(B11[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B12[b,r])))) + dot(B11[b,r],dmul(Gamma{11},dmul(Gamma{5}, B12[b,r]))) + dot(B11[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B11[b,r]))))
+dot(B12[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B11[b,r])))) + dot(B12[b,r],dmul(Gamma{11},dmul(Gamma{5},-B11[b,r]))) + dot(B12[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B12[b,r]))))
+dot(B21[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B22[b,r])))) + dot(B21[b,r],dmul(Gamma{11},dmul(Gamma{5}, B22[b,r]))) + dot(B21[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B21[b,r]))))
+dot(B22[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B21[b,r])))) + dot(B22[b,r],dmul(Gamma{11},dmul(Gamma{5},-B21[b,r]))) + dot(B22[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B22[b,r]))))
+dot(B31[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B32[b,r])))) + dot(B31[b,r],dmul(Gamma{11},dmul(Gamma{5}, B32[b,r]))) + dot(B31[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B31[b,r]))))
+dot(B32[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B31[b,r])))) + dot(B32[b,r],dmul(Gamma{11},dmul(Gamma{5},-B31[b,r]))) + dot(B32[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B32[b,r])))) )/(3*prod(lp.iL[1:3]))
end end end
end
return cf
end
f1_fun() = (norm2(C11) + norm2(C21) + norm2(C31) + norm2(C12) + norm2(C22) + norm2(C32))
k1_fun() = (dot(C11,dmul(Gamma{6},imm(C12))) + dot(C11,dmul(Gamma{7},-C12)) + dot(C11,dmul(Gamma{8}, imm(C11)))
+dot(C12,dmul(Gamma{6},imm(C11))) + dot(C12,dmul(Gamma{7}, C11)) + dot(C12,dmul(Gamma{8},mimm(C12)))
+dot(C21,dmul(Gamma{6},imm(C22))) + dot(C21,dmul(Gamma{7},-C22)) + dot(C21,dmul(Gamma{8}, imm(C21)))
+dot(C22,dmul(Gamma{6},imm(C21))) + dot(C22,dmul(Gamma{7}, C21)) + dot(C22,dmul(Gamma{8},mimm(C22)))
+dot(C31,dmul(Gamma{6},imm(C32))) + dot(C31,dmul(Gamma{7},-C32)) + dot(C31,dmul(Gamma{8}, imm(C31)))
+dot(C32,dmul(Gamma{6},imm(C31))) + dot(C32,dmul(Gamma{7}, C31)) + dot(C32,dmul(Gamma{8},mimm(C32))) )/3
function compute_correlators()
println("Computing correlators")
println("Computing f1...")
file = open("./output/"*basename(parsed_args["c"])*".f1","w+")
global f1 = f1_fun()
println(file," ",f1)
flush(file)
close(file)
println("Computing k1...")
file = open("./output/"*basename(parsed_args["c"])*".k1","w+")
global k1 = k1_fun()
println(file," ",k1)
flush(file)
close(file)
println("Computing fP...")
file = open("./output/"*basename(parsed_args["c"])*".fP","w+")
global fP = fP_fun()
println(file," ",fP)
flush(file)
close(file)
println("Computing fA...")
file = open("./output/"*basename(parsed_args["c"])*".fA","w+")
global fA = fA_fun()
flush(file)
close(file)
println("Computing gP...")
file = open("./output/"*basename(parsed_args["c"])*".gP","w+")
global gP = gp_fun()
flush(file)
close(file)
println("Computing gA...")
file = open("./output/"*basename(parsed_args["c"])*".gA","w+")
global gA = ga_fun()
flush(file)
close(file)
println("Computing kV...")
file = open("./output/"*basename(parsed_args["c"])*".kV","w+")
global kv = kv_fun()
flush(file)
close(file)
println("Computing lV...")
file = open("./output/"*basename(parsed_args["c"])*".lV","w+")
global lv = lv_fun()
flush(file)
close(file)
println("Computing kT...")
file = open("./output/"*basename(parsed_args["c"])*".kT","w+")
global kt = kt_fun()
flush(file)
close(file)
println("Computing lT...")
file = open("./output/"*basename(parsed_args["c"])*".lT","w+")
global lt = lt_fun()
flush(file)
close(file)
end

160
src/io.jl Normal file
View file

@ -0,0 +1,160 @@
using LatticeGPU
using TOML
using TimerOutputs
using ArgParse
using CUDA
"""
function read_input()
Stores as global variables 'parsed_args' (info from the command line) and 'params' (info from the input file)
"""
function read_input()
global parsed_args = parse_commandline()
println("--------------------------------------------------")
println("Reading input file from:", parsed_args["i"], "...")
global params = TOML.parsefile(parsed_args["i"])
return nothing
end
function parse_commandline()
s = ArgParseSettings()
@add_arg_table s begin
"-i"
help = "Input parameters file"
required = true
arg_type = String
"-c"
help = "Gauge configuration file"
required = true
arg_type = String
"--cern"
help = "Config written with the export_cnfg_cern() convention"
action = :store_true
end
return parse_args(s)
end
"""
function load_gauge_field()
Returns the gauge field and computes the Csw term
"""
function load_gauge_field()
println("Reading gauge field from: ", parsed_args["c"], "...")
if !parsed_args["cern"]
U,_ = read_cnfg(parsed_args["c"])
else
U = read_cnfg_cern(parsed_args["c"],lp)
end
Csw!(dws, U, gp, lp)
return U
end
function read_cnfg_cern(path::String,lp::SpaceParm)
U = vector_field(SU3{Float64}, lp);
file = open(path)
for t in 1:lp.iL[4]
for i in 1:lp.iL[1]
for j in 1:lp.iL[2]
for k in 1:lp.iL[3]
for d in [4,1,2,3]
f,r = point_index(CartesianIndex((i,j,k,t)),lp)
#a11 !!
re11 = read(file,Float64)
co11 = read(file,Float64)
#a12 !!
re12 = read(file,Float64)
co12 = read(file,Float64)
#a13 !!
re13 = read(file,Float64)
co13 = read(file,Float64)
#a21 !!
re21 = read(file,Float64)
co21 = read(file,Float64)
#a22 !!
re22 = read(file,Float64)
co22 = read(file,Float64)
#a23 !!
re23 = read(file,Float64)
co23 = read(file,Float64)
#a31
re31 = read(file,Float64)
co31 = read(file,Float64)
#a32
re32 = read(file,Float64)
co32 = read(file,Float64)
#a33
re33 = read(file,Float64)
co33 = read(file,Float64)
CUDA.@allowscalar (U[f,d,r] = SU3{Float64}(re11 + im*co11, re12 + im*co12, re13 + im*co13,
re21 + im*co21, re22 + im*co22, re23 + im*co23))
end
end
end
end
end
length(read(file)) == (prod(lp.iL[1:3])*4*8*9*2) ? nothing : error("File not fully read")
close(file)
return U
end
"""
function load_structs()
Stores in global variables the needed structures, i.e. lp, gp, dpar, dws, ymws
"""
function load_structs()
global lp = SpaceParm{4}(tuple(params["Space"]["size"]...), tuple(params["Space"]["blocks"]...),BC_SF_ORBI, (0,0,0,0,0,0))
global gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(params["Space"]["cG"],0.0),params["Space"]["phiT"],lp.iL);
global dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,params["Fermion"]["csw"],ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,params["Fermion"]["ct"]);
global dws = DiracWorkspace(SU3fund,Float64,lp);
global ymws = YMworkspace(SU3,Float64,lp);
println("Parameters:")
println("Lattice size: ", lp.iL)
println("Phi0 = ", params["Space"]["phi0"])
println("PhiT = ", params["Space"]["phiT"])
println("cG = ", gp.cG[1])
println("kappa = ", params["Fermion"]["kappa"])
println("theta = ", params["Fermion"]["theta"])
println("csw = ", dpar.csw)
println("ct = ", dpar.ct)
println("tolerance = ", params["Solver"]["tolerance"])
println("maxiter = ", params["Solver"]["maxiter"])
return nothing
end

77
src/props.jl Normal file
View file

@ -0,0 +1,77 @@
using LatticeGPU
using TOML
using TimerOutputs
using ArgParse
using CUDA
"""
function compute_propagators()
Computes the propagators for the first two components of spin.
The propagators are stored in Acs (t=0 to bulk), Bcs(bulk to t=T) and Ccs(t=0 to t=T)
"""
function compute_propagators()
println("--------------------------------------------------")
psi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
println("Computing propagators from t=0 to the bulk")
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A11 = Array(psi)
global C11 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
println("CG converged for c=2 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A21 = Array(psi)
global C21 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 1)
println("CG converged for c=3 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A31 = Array(psi)
global C31 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 2)
println("CG converged for c=1 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A12 = Array(psi)
global C12 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 2)
println("CG converged for c=2 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A22 = Array(psi)
global C22 = bndtobnd(psi, U, dpar, dws, lp)
niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 2)
println("CG converged for c=3 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global A32 = Array(psi)
global C32 = bndtobnd(psi, U, dpar, dws, lp)
println("Computing propagators from t=T to the bulk")
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B11 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
println("CG converged for c=2 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B21 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 1)
println("CG converged for c=3 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B31 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 2)
println("CG converged for c=1 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B12 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 2)
println("CG converged for c=2 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B22 = Array(psi)
niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 2)
println("CG converged for c=3 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
global B32 = Array(psi)
println("--------------------------------------------------")
return nothing
end