Runnable version

input/FerFlow.in

 [Run]
 user    = "Lattice_user"
-name    = "Template_run"
+name    = "Run1L32"
 
 [Space]
-size       = [8,8,8,8]
+size       = [32,32,32,32]
 blocks     = [4,4,4,4]
 
 [HMC]
@@ -13,10 +13,9 @@ dmeas = 5
 eps = 0.01
 ns = 30
 
-
 [Fermion]
-beta   = 8.4044000000000008e+00
-kappa  = 1.3247670000000000e-01
+beta   = 6.0
+kappa  = 0.1348
 theta  = 0.0
 theta_t= 3.1415926535897
 csw    = 1.76923076923077
@@ -28,4 +27,4 @@ N_noise = 12
 Flow_times = [2.47,3.86]
 Nsaves = 25
 tolerance = 1.0e-13
-maxiter   = 5000
+maxiter   = 50000

main.jl

 using Pkg
-Pkg.activate("./");
+Pkg.activate("/home/fperez/Git/LGPU_fork_ferflow")
 
 using LatticeGPU
 using TOML
@@ -10,10 +10,10 @@ using InteractiveUtils
 using BDIO
 using MD5
 
+reset_timer!()
+
 include("./src/io.jl")
 include("./src/meas.jl")
-include("./src/corr.jl")
-
 
 @timeit "Input reading and space allocation" begin
     read_input()
@@ -23,7 +23,7 @@ include("./src/corr.jl")
     load_fields()
 end
 
-thermalize()
+@timeit "Thermalization" thermalize()
 
 while isfile("run_on")
 

output/Run1L32.log

+Running Fermass.jl by Lattice_user. Name of the run: Run1L32
+
+Calling: main.jl -i ./input/FerFlow.in 
+
+Version info:
+Julia Version 1.8.5
+Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
+Platform Info:
+  OS: Linux (x86_64-linux-gnu)
+  CPU: 128 × AMD EPYC 7452 32-Core Processor
+  WORD_SIZE: 64
+  LIBM: libopenlibm
+  LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
+  Threads: 1 on 128 virtual cores
+
+Reading input file from:./input/FerFlow.in...
+
+Parameters:
+Lattice size:  (32, 32, 32, 32)
+kappa        = 0.1348
+theta        = 0.0
+csw          = 1.76923076923077
+tolerance    = 1.0e-13
+maxiter      = 50000
+Flow times   = [2.47, 3.86]
+N noise      = 12
+
+ Starting simulation
+
+Thermalizing...
+Performing 5 HMC steps...
+Updated to Gauge configuration n1
+
+Measuring 2pt...
+Finished measuring 2pt
+
+Measuring 1pt at t = 2.47...
+Measuring 1pt at t = 3.86...
+Finished measuring 1pt
+
+Creating new BDIO output file ./output/Run1L32.bdio
+Performing 5 HMC steps...
+Updated to Gauge configuration n2
+
+Measuring 2pt...
+Finished measuring 2pt
+
+Measuring 1pt at t = 2.47...
+Measuring 1pt at t = 3.86...
+Finished measuring 1pt
+
+
+Appending output to ./output/Run1L32.bdio
+
+Performing 5 HMC steps...
+Updated to Gauge configuration n3
+
+Measuring 2pt...
+Finished measuring 2pt
+
+Measuring 1pt at t = 2.47...
+Measuring 1pt at t = 3.86...
+Finished measuring 1pt
+
+
+Appending output to ./output/Run1L32.bdio
+
+Performing 5 HMC steps...
+Updated to Gauge configuration n4
+
+Measuring 2pt...
+Finished measuring 2pt
+
+Measuring 1pt at t = 2.47...
+Measuring 1pt at t = 3.86...
+Finished measuring 1pt
+
+
+Appending output to ./output/Run1L32.bdio
+
+Performing 5 HMC steps...
+Updated to Gauge configuration n5
+
+Measuring 2pt...

src/io.jl

@@ -50,9 +50,9 @@ 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 lp = SpaceParm{4}(tuple(params["Space"]["size"]...), tuple(params["Space"]["blocks"]...),0, (0,0,0,0,0,0))
+    global gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(0.0,0.0),(0.0,0.0),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),1.0);
     global dws = DiracWorkspace(SU3fund,Float64,lp);
     global ymws = YMworkspace(SU3,Float64,lp);
     global int = wfl_rk3(Float64, 0.1, 1.0)
@@ -82,13 +82,9 @@ function write_log()
 
     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)
@@ -101,7 +97,7 @@ end
 
 function save_data()
 
-    ihdr = [convert(Int32,2708686513)]
+    ihdr = [convert(Int32,0608686513)]
     fname = "./output/"*params["Run"]["name"]*".bdio"
 
     if isfile(fname)

src/meas.jl

@@ -6,12 +6,14 @@ using CUDA
 
 function load_fields()
     global U = vector_field(SU3{Float64}, lp)
+    fill!(U,one(SU3{Float64}));
+
     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_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));
@@ -19,22 +21,28 @@ function load_fields()
     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));
+    global Quark_cond = Array{Complex{Float64}}(undef,N_noise,length(flow_times));
+    return nothing
 end
 
 
-function themalize()
-    println(log_file,"Thermalizing")
+function thermalize()
+
+    println(log_file,"\n Starting simulation\n")
+	
+    println(log_file,"Thermalizing...")
     flush(log_file)
 
     for _ in 1:params["HMC"]["nth"]
         HMC!(U,intsch,lp,gp,ymws)
     end
+
+    return nothing
 end
 
 function Gauge_update()
 
-    println(log_file,"Performing "*string(params["HMC"]["dmeas"])*" HMC steps")
+    println(log_file,"Performing "*string(params["HMC"]["dmeas"])*" HMC steps...")
 
     flush(log_file)
 
@@ -43,7 +51,9 @@ function Gauge_update()
         acc ? nothing : println(log_file,"Rejected HCM step")
     end
 
-    MCid += 1
+
+    global MCid = MCid + 1
+
 
     println(log_file,"Updated to Gauge configuration n"*string(MCid))
 
@@ -62,23 +72,23 @@ Stores the two point function in the variables 'pp_corr_t0', 'ap_corr_t0', 'pp_c
 """
 function two_pt()
 
+        println(log_file,"\nMeasuring 2pt...")
+        flush(log_file)
     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)))
+        ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
 
-        pp_corr_t0 .= zero(Float64)
-        ap_corr_t0 .= zero(ComplexF64)
+        pp_corr_t0[:,noi] .= zero(Float64)
+        ap_corr_t0[:,noi] .= 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
+                pp_corr_t0[t,noi] += pp_density[b,r]
+                ap_corr_t0[t,noi] += ap_density[b,r]
             end end end
         end
 
@@ -87,20 +97,19 @@ function two_pt()
         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)
+            ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
+            pp_corr_t[:,noi,fl] .= zero(Float64)
+            ap_corr_t[:,noi,fl] .= 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
+                    pp_corr_t[t,noi,fl] += pp_density[b,r]
+                    ap_corr_t[t,noi,fl] += ap_density[b,r]
                 end end end
             end
         end
@@ -120,12 +129,12 @@ 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)*"...")
+            println(log_file,"Measuring 1pt at t = "*string(flow_times[fl])*"...")
             flush(log_file)
 
+    	for noi in 1:N_noise
+
             pfrandomize!(psi,lp)
 
             backflow(psi, U, flow_times[fl], params["Measurements"]["Nsaves"], gp, dpar, lp, ymws, dws)
@@ -138,7 +147,7 @@ function one_pt()
             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))
+            Quark_cond[noi,fl] = sum(dot.(dws.sp,psi) - params["Fermion"]["Cfl"] .* norm2.(dws.sp))
 
             copyto!(U,U_CPU)