Initial commit

README.md

+

input/FerFlow.in

+[Run]
+user    = "Lattice_user"
+name    = "Template_run"
+
+[Space]
+size       = [8,8,8,8]
+blocks     = [4,4,4,4]
+
+[HMC]
+nth = 5000
+beta = 6.0
+dmeas = 5
+eps = 0.01
+ns = 30
+
+
+[Fermion]
+beta   = 8.4044000000000008e+00
+kappa  = 1.3247670000000000e-01
+theta  = 0.0
+theta_t= 3.1415926535897
+csw    = 1.76923076923077
+ct     = 9.8718528389891003e-01
+Cfl    = 0.5734
+
+[Measurements]
+N_noise = 12
+Flow_times = [2.47,3.86]
+Nsaves = 25
+tolerance = 1.0e-13
+maxiter   = 5000

main.jl

+using Pkg
+Pkg.activate("./");
+
+using LatticeGPU
+using TOML
+using TimerOutputs
+using ArgParse
+using CUDA
+using InteractiveUtils
+using BDIO
+using MD5
+
+include("./src/io.jl")
+include("./src/meas.jl")
+include("./src/corr.jl")
+
+
+@timeit "Input reading and space allocation" begin
+    read_input()
+    load_structs()
+    global log_file = open("./output/"*params["Run"]["name"]*".log","w+")
+    write_log()
+    load_fields()
+end
+
+thermalize()
+
+while isfile("run_on")
+
+    @timeit "HMC" Gauge_update()
+
+    @timeit "Two point function" two_pt()
+
+    @timeit "One point function" one_pt()
+
+    @timeit "Saving" save_data()
+
+end
+
+print_timer(log_file)
+
+flush(log_file)
+close(log_file)

src/io.jl

+using LatticeGPU
+using TOML
+using ArgParse
+using CUDA
+using BDIO
+using InteractiveUtils
+using MD5
+
+
+"""
+    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()
+    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"
+        arg_type = String
+
+        "--cern"
+        help = "Config written with the export_cnfg_cern() convention"
+        action = :store_true
+
+    end
+
+    return parse_args(s)
+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]) + ((i==4)*im*params["Fermion"]["theta_t"]/lp.iL[i])),4),0.0,params["Fermion"]["ct"]);
+    global dws = DiracWorkspace(SU3fund,Float64,lp);
+    global ymws = YMworkspace(SU3,Float64,lp);
+    global int = wfl_rk3(Float64, 0.1, 1.0)
+    global intsch = omf4(Float64,params["HMC"]["eps"], params["HMC"]["ns"]);
+
+    global flow_times = params["Measurements"]["Flow_times"]
+    global N_noise = params["Measurements"]["N_noise"]
+    global MCid = 0
+
+    return nothing
+end
+
+function write_log()
+
+    println(log_file,"Running Fermass.jl by ", params["Run"]["user"],". Name of the run: ",params["Run"]["name"])
+    println(log_file,"")
+
+    print(log_file,"Calling: ")
+    print(log_file,PROGRAM_FILE*" "); for x in ARGS; print(log_file,x*" "); end
+    println(log_file,"\n")
+
+    println(log_file,"Version info:")
+    versioninfo(log_file)
+    println(log_file,"")
+
+    println(log_file,"Reading input file from:", parsed_args["i"], "...\n")
+
+    println(log_file,"Parameters:")
+    println(log_file,"Lattice size:  ", lp.iL)
+    println(log_file,"Phi0         = ", params["Space"]["phi0"])
+    println(log_file,"PhiT         = ", params["Space"]["phiT"])
+    println(log_file,"cG           = ", gp.cG[1])
+    println(log_file,"kappa        = ", params["Fermion"]["kappa"])
+    println(log_file,"theta        = ", params["Fermion"]["theta"])
+    println(log_file,"csw          = ", dpar.csw)
+    println(log_file,"ct           = ", dpar.ct)
+    println(log_file,"tolerance    = ", params["Measurements"]["tolerance"])
+    println(log_file,"maxiter      = ", params["Measurements"]["maxiter"])
+    println(log_file,"Flow times   = ", flow_times)
+    println(log_file,"N noise      = ", N_noise)
+
+    flush(log_file)
+    return nothing
+end
+
+
+function save_data()
+
+    ihdr = [convert(Int32,2708686513)]
+    fname = "./output/"*params["Run"]["name"]*".bdio"
+
+    if isfile(fname)
+        fb = BDIO_open(fname, "a")
+        println(log_file,"\nAppending output to "*fname*"\n")
+    else
+        fb = BDIO_open(fname, "w", "BDIO output from sfcf.jl")
+        println(log_file,"Creating new BDIO output file "*fname)
+
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 14)
+        BDIO_write!(fb, ihdr)
+        BDIO_write_hash!(fb)
+
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 1)
+        BDIO_write!(fb, [convert(Int32, 4)])
+        BDIO_write!(fb, [convert(Int32, lp.iL[i]) for i in 1:4])
+        BDIO_write!(fb, [dpar.m0, dpar.csw, dpar.ct,params["Fermion"]["Cfl"]])
+        BDIO_write!(fb, [dpar.th[i] for i in 1:4])
+        BDIO_write!(fb, flow_times)
+        BDIO_write!(fb, [N_noise,])
+
+        BDIO_write_hash!(fb)
+    end
+
+
+    BDIO_start_record!(fb, BDIO_BIN_GENERIC, 8)
+
+    for fl in 1:length(flow_times)
+        for noi in 1:N_noise
+            BDIO_write!(fb,pp_corr_t0)
+            BDIO_write!(fb,ap_corr_t0)
+
+            BDIO_write!(fb,pp_corr_t)
+            BDIO_write!(fb,ap_corr_t)
+            BDIO_write!(fb,Quark_cond)
+        end
+    end
+
+    BDIO_write_hash!(fb)
+
+    BDIO_close!(fb)
+
+    return nothing
+end
+
+
+
+
+"""
+    function load_gauge_field()
+
+Returns the gauge field and computes the Csw term
+"""
+function load_gauge_field()
+    if !parsed_args["cern"]
+    println(log_file,"\nReading gauge field with LGPU native convention from: ", parsed_args["c"], "...")
+        U,_ = read_cnfg(parsed_args["c"])
+    else
+    println(log_file,"\nReading gauge field with CERN convention from: ", parsed_args["c"], "...")
+        U = read_cnfg_cern(parsed_args["c"],lp)
+    end
+
+    f = open(parsed_args["c"],"r")
+    println(log_file,"MD5 checksum of gauge config: ",bytes2hex(md5(f)),"\n")
+    close(f)
+
+    Csw!(dws, U, gp, lp)
+    return U
+end
+
+function read_cnfg_cern(path::String,lp::SpaceParm)
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    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)
+
+                        (Ucpu[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 CuArray(Ucpu)
+end

src/meas.jl

+using LatticeGPU
+using TOML
+using TimerOutputs
+using ArgParse
+using CUDA
+
+function load_fields()
+    global U = vector_field(SU3{Float64}, lp)
+    global U_CPU = Array(U)
+    global psi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+    global psi_CPU = Array(psi)
+
+    global pp_density = Array{Float64}(undef,lp.bsz,lp.rsz));
+    global ap_density = Array{Complex{Float64}}(undef,lp.bsz,lp.rsz));
+
+    global pp_corr_t0 = fill(zero(Float64),(lp.iL[4],N_noise));
+    global ap_corr_t0 = fill(zero(ComplexF64),(lp.iL[4],N_noise));
+
+    global pp_corr_t = fill(zero(Float64),(lp.iL[4],N_noise,length(flow_times)));
+    global ap_corr_t = fill(zero(ComplexF64),(lp.iL[4],N_noise,length(flow_times)));
+
+    global Quark_cond = Array{Complex{Float64}}(undef,noise,length(flow_times));
+end
+
+
+function themalize()
+    println(log_file,"Thermalizing")
+    flush(log_file)
+
+    for _ in 1:params["HMC"]["nth"]
+        HMC!(U,intsch,lp,gp,ymws)
+   end
+end
+
+function Gauge_update()
+
+    println(log_file,"Performing "*string(params["HMC"]["dmeas"])*" HMC steps")
+
+    flush(log_file)
+
+    for _ in 1:params["HMC"]["dmeas"]
+        _,acc = HMC!(U,intsch,lp,gp,ymws)
+        acc ? nothing : println(log_file,"Rejected HCM step")
+    end
+
+    MCid += 1
+
+    println(log_file,"Updated to Gauge configuration n"*string(MCid))
+
+    Csw!(dws,U,gp,lp)
+    U_CPU .= Array(U)
+
+    flush(log_file)
+    return nothing
+
+end
+
+"""
+function two_pt()
+
+Stores the two point function in the variables 'pp_corr_t0', 'ap_corr_t0', 'pp_corr_t' and 'ap_corr_t'
+"""
+function two_pt()
+
+    for noi in 1:N_noise
+
+        println(log_file,"Measuring 2pt noise source "*string(noi)*"...")
+        flush(log_file)
+
+        propagator!(psi,U, dpar, dws, lp,params["Measurements"]["maxiter"], params["Measurements"]["tolerance"],1)
+
+        pp_density .= Array(norm2.(psi))
+        ap_density .= Array(dot.(psi,dmul(Gamma{4},psi)))
+
+        pp_corr_t0 .= zero(Float64)
+        ap_corr_t0 .= zero(ComplexF64)
+        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]
+                b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+                pp_corr_t0[t,noi] += pp_density
+                ap_corr_t0[t,noi] += ap_density
+            end end end
+        end
+
+
+
+        delta_t = 0.0
+        neps = int.eps_ini
+        for fl in 1:length(flow_times)
+
+            delta_t = flow_times[fl] - delta_t
+            _,epslist = flw_adapt(U, psi, int, delta_t, neps, gp, dpar, lp, ymws, dws)
+            neps = epslist[end]
+
+            pp_density .= Array(norm2.(psi))
+            ap_density .= Array(dot.(psi,dmul(Gamma{4},psi)))
+            pp_corr_t .= zero(Float64)
+            ap_corr_t .= zero(ComplexF64)
+            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]
+                    b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+                    pp_corr_t[t,noi,fl] += pp_density
+                    ap_corr_t[t,noi,fl] += ap_density
+                end end end
+            end
+        end
+
+        copyto!(U,U_CPU)
+    end
+
+    println(log_file,"Finished measuring 2pt\n")
+    flush(log_file)
+
+    return nothing
+end
+
+
+
+function one_pt()
+
+    copyto!(U,U_CPU)
+
+    for noi in 1:N_noise
+        for fl in 1:length(flow_times)
+
+            println(log_file,"Measuring 1pt noise source "*string(noi)*"...")
+            flush(log_file)
+
+            pfrandomize!(psi,lp)
+
+            backflow(psi, U, flow_times[fl], params["Measurements"]["Nsaves"], gp, dpar, lp, ymws, dws)
+
+            psi_CPU .= Array(psi)
+
+            dws.sp .= dmul.(Gamma{5},psi)
+            g5Dw!(psi,U,dws.sp,dpar,dws,lp)
+
+            CG!(psi,U,DwdagDw!,dpar,lp,dws,params["Measurements"]["maxiter"], params["Measurements"]["tolerance"])
+
+            copyto!(dws.sp,psi_CPU)
+            Quark_cond[noi,fl] = sum(dot.(dws.sp,psi) - params["Fermion"]["Cfl"].*norm2.(dws.sp))
+
+            copyto!(U,U_CPU)
+
+        end
+    end
+    println(log_file,"Finished measuring 1pt\n")
+
+    return nothing
+end
+
+function gauge_meas()
+
+    println(log_file,"\nGauge meassurements:\n")
+    println(log_file,"Plaquette: ", plaquette(U, lp, gp, ymws))
+    println(log_file,"Topological charge: ", Qtop(U, gp, lp, ymws))
+    sfc = sfcoupling(U, lp, gp, ymws)
+    flush(log_file)
+
+    return nothing
+end