From 3cb578d0331403d0c8738741c47619369719b681 Mon Sep 17 00:00:00 2001 From: "Fernando P. Panadero" Date: Wed, 21 Feb 2024 17:05:34 +0100 Subject: [PATCH] Split of the code in src --- sfcf.jl | 358 ++------------------------------------------------- src/corr.jl | 213 ++++++++++++++++++++++++++++++ src/io.jl | 160 +++++++++++++++++++++++ src/props.jl | 77 +++++++++++ 4 files changed, 459 insertions(+), 349 deletions(-) create mode 100644 src/corr.jl create mode 100644 src/io.jl create mode 100644 src/props.jl diff --git a/sfcf.jl b/sfcf.jl index bead746..39b6c27 100644 --- a/sfcf.jl +++ b/sfcf.jl @@ -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") diff --git a/src/corr.jl b/src/corr.jl new file mode 100644 index 0000000..3a87321 --- /dev/null +++ b/src/corr.jl @@ -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 diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 0000000..79da828 --- /dev/null +++ b/src/io.jl @@ -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 diff --git a/src/props.jl b/src/props.jl new file mode 100644 index 0000000..a998c71 --- /dev/null +++ b/src/props.jl @@ -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