From 713812c3b90756d7619c385d6f938cc36ef3f999 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fernando=20P=C3=A9rez=20Panadero?= Date: Mon, 19 Feb 2024 17:20:28 +0100 Subject: [PATCH] Main script and dirs. Only four correlation functions. --- input/sfcf.in | 21 +++++ sfcf.jl | 211 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 232 insertions(+) create mode 100644 input/sfcf.in create mode 100644 sfcf.jl diff --git a/input/sfcf.in b/input/sfcf.in new file mode 100644 index 0000000..8e54bd6 --- /dev/null +++ b/input/sfcf.in @@ -0,0 +1,21 @@ +[Run] +user = "fperez" +name = "SFtest_Setup10" + +[Space] +size = [8,8,8,8] +blocks = [4,4,4,4] +phi0 = [-1.047197551196598, 0.0] # Boundary fields for SF +phiT = [-3.141592653589793, 1.047197551196598] +cG = 1.0 # Correction for boundaries (SF only) + +[Fermion] +beta = 6.0 #Used only if cG/csw/ct are not input +kappa = 0.1 +theta = 0.5 +csw = 1.0 +ct = 0.9 + +[Solver] +tolerance = 1.0e-10 +maxiter = 1000 diff --git a/sfcf.jl b/sfcf.jl new file mode 100644 index 0000000..0e12d45 --- /dev/null +++ b/sfcf.jl @@ -0,0 +1,211 @@ +using LatticeGPU +using TOML +using TimerOutputs +using ArgParse +using CUDA + +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 + + 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"]) + +cG = try + params["Space"]["cG"] +catch + try + println("Using default value for cG") + 1.0 - 0.089*6/params["Fermion"]["beta"] + catch + error("Either cG of beta has to be defined in the input file") + end +end + +ct = try + params["Fermion"]["ct"] +catch + try + error("ct not defined") + println("Using default value for ct") + catch + error("Either ct of beta has to be defined in the input file") + end +end + +csw = try + params["Fermion"]["csw"] +catch + try + error("csw not defined") + println("Using default value for csw") + catch + error("Either csw of beta has to be defined in the input file") + end +end + +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,(cG,0.0),params["Space"]["phiT"],lp.iL); +dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,csw,ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,ct); + +dws = DiracWorkspace(SU3fund,Float64,lp); +ymws = YMworkspace(SU3,Float64,lp); + +println("Reading gauge field from: ", parsed_args["c"], "...") +U,_ = read_cnfg(parsed_args["c"]) +Csw!(dws, U, gp, lp) +println("--------------------------------------------------") + +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) + +f1 = 0.0 + +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 s=1 and c=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr)) +A11 = Array(psi) +f1 += norm2(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) +f1 += norm2(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) +f1 += norm2(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) +f1 += norm2(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) +f1 += norm2(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) +f1 += norm2(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 s=1 and c=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 fP...") +file = open(basename(parsed_args["i"])*".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]) + end end end + println(file,t," ",fP[t]) +end +flush(file) +close(file) + +println("Computing fA...") +file = open(basename(parsed_args["i"])*".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],A13[b,r]) + dot(A21[b,r],A23[b,r]) + dot(A31[b,r],A33[b,r]) + dot(A12[b,r],A14[b,r]) + dot(A22[b,r],A24[b,r]) + dot(A32[b,r],A34[b,r]) + end end end + println(file,t," ",fA[t]) +end +flush(file) +close(file) + +println("Computing gP...") +file = open(basename(parsed_args["i"])*".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]) + println(file,t," ",gP[t]) + end end end +end +flush(file) +close(file) + + +println("Computing gA...") +file = open(basename(parsed_args["i"])*".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],B13[b,r]) + dot(B21[b,r],B23[b,r]) + dot(B31[b,r],B33[b,r]) + dot(B12[b,r],B14[b,r]) + dot(B22[b,r],B24[b,r]) + dot(B32[b,r],B34[b,r]) + end end end + println(file,t," ",gA[t]) +end +flush(file) +close(file)