diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 5cf7063..96c9c0c 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -84,20 +84,23 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, cG, ztw, lp::SpaceParm{N,M,B @inbounds begin b = Int64(CUDA.threadIdx().x) r = Int64(CUDA.blockIdx().x) + TWP = (b==1) && (r==1) + it = point_time((b, r), lp) - + IBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && + ( (it == 1) || (it == lp.iL[end])) ) + S = zero(eltype(plx)) ipl = 0 for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) - SFBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && - ( (it == 1) || (it == lp.iL[end])) ) && (id1 == N) + SFBND = IBND && (id1 == N) for id2 = 1:id1-1 bu2, ru2 = up((b, r), id2, lp) ipl = ipl + 1 - if SFBND && (it == lp.iL[end]) + if SFBND && (it == lp.iL[end]) gt1 = Ubnd[id2] else gt1 = U[bu1,id2,ru1] @@ -106,7 +109,11 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, cG, ztw, lp::SpaceParm{N,M,B if SFBND S += cG*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) else - S += ztw[ipl]*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) + if TWP + S += ztw[ipl]*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) + else + S += tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) + end end end end @@ -122,6 +129,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw, b = Int64(CUDA.threadIdx().x) r = Int64(CUDA.blockIdx().x) + TWP = (b==1) && (r==1) it = point_time((b,r), lp) @inbounds begin @@ -153,12 +161,17 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw, frc1[b ,id2, r ] += X frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) else - X = projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) - + if TWP + X = projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) + frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2) + frc2[bu2,id1,ru2] += projalg(ztw,g2*g1) + else + X = projalg(U[b,id1,r]*g1/U[b,id2,r]) + frc2[bu1,id2,ru1] -= projalg(g1*g2) + frc2[bu2,id1,ru2] += projalg(g2*g1) + end frc1[b ,id1, r ] -= X frc1[b ,id2, r ] += X - frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2) - frc2[bu2,id1,ru2] += projalg(ztw,g2*g1) end end diff --git a/src/main/times.jl b/src/main/times.jl index 38d9a18..9286413 100644 --- a/src/main/times.jl +++ b/src/main/times.jl @@ -7,7 +7,7 @@ Pkg.activate("../..") using LatticeGPU # Set lattice/block size -ntwist = (0,0,0,0,0,0) +ntwist = (1,0,0,0,0,0) lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist) println("Space Parameters: ", lp) @@ -41,10 +41,10 @@ println("Time to take the configuration to memory: ") # Set gauge parameters # FIRST SET: Wilson action/flow println("\n## WILSON ACTION/FLOW TIMES") -gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.66666666) +gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0) println("Gauge Parameters: ", gp) -flwint = zfl_rk3(PREC, 0.005, 1.0E-6) +flwint = wfl_rk3(PREC, 0.005, 1.0E-6) println(flwint) println("Initial Action: ") @@ -54,7 +54,7 @@ println("Initial Action: ") HMC!(U,int.eps,1,lp, gp, ymws, noacc=true) pl = Vector{Float64}() -for i in 1:100 +for i in 1:4 @time dh, acc = HMC!(U,int,lp, gp, ymws) println("# HMC: ", acc, " ", dh) push!(pl, plaquette(U,lp, gp, ymws)) @@ -69,7 +69,7 @@ println("Action: ", gauge_action(U, lp, gp, ymws)) flw(U, flwint, 1, gp, lp, ymws) println("Time for 100 steps of RK3 flow integrator: ") -@time flw(U, flwint, 400, gp, lp, ymws) +@time flw(U, flwint, 100, gp, lp, ymws) eoft = Eoft_plaq(U, gp, lp, ymws) eoft = Eoft_clover(U, gp, lp, ymws) qtop = Qtop(U, gp, lp, ymws)