lattice/latticegpu.jl
1
0
Fork 0
mirror of https://igit.ific.uv.es/alramos/latticegpu.jl.git synced 2025-05-14 19:23:42 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Full support for twisted bc

Browse source
This commit is contained in:
Alberto Ramos 2021-10-31 08:37:53 +01:00
parent 26b5dfab01
commit 410a6d7b25
9 changed files with 214 additions and 100 deletions
Download patch file Download diff file Expand all files Collapse all files

24
src/Groups/AlgebraSU3.jl
View file

@ -24,6 +24,30 @@ function projalg(a::SU3{T}) where T <: AbstractFloat
sr3ov2*(ditr)) sr3ov2*(ditr))
end end
function projalg(z::Complex{T}, a::SU3{T}) where T <: AbstractFloat
sr3ov2::T = 0.866025403784438646763723170752
zu11 = z*a.u11
zu12 = z*a.u12
zu13 = z*a.u13
zu21 = z*a.u21
zu22 = z*a.u22
zu23 = z*a.u23
ditr = ( imag(zu11) + imag(zu22) - 2.0*imag(z*conj(a.u11*a.u22 - a.u12*a.u21)) )/3.0
m12 = (zu12 - conj(zu21))/2.0
m13 = (zu13 - conj(z)*(a.u12*a.u23 - a.u13*a.u22) )/2.0
m23 = (zu23 - conj(z)*(a.u13*a.u21 - a.u11*a.u23) )/2.0
return SU3alg{T}(imag( m12 ), imag( m13 ), imag( m23 ),
real( m12 ), real( m13 ), real( m23 ),
(imag(zu11)-imag(zu22))/2.0,
sr3ov2*(ditr))
end
projalg(z::Int, a::SU3{T}) where T <: AbstractFloat = projalg(a::SU3{T})
dot(a::SU3alg{T},b::SU3alg{T}) where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 + a.t4*b.t4 + a.t5*b.t5 + a.t6*b.t6 + a.t7*b.t7 + a.t8*b.t8 dot(a::SU3alg{T},b::SU3alg{T}) where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 + a.t4*b.t4 + a.t5*b.t5 + a.t6*b.t6 + a.t7*b.t7 + a.t8*b.t8
norm2(a::SU3alg{T}) where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2 norm2(a::SU3alg{T}) where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2
norm(a::SU3alg{T}) where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2) norm(a::SU3alg{T}) where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2)

2
src/Groups/M3x3.jl
View file

@ -117,7 +117,7 @@ function projalg(a::M3x3{T}) where T <: AbstractFloat
sr3ov2::T = 0.866025403784438646763723170752 sr3ov2::T = 0.866025403784438646763723170752
ditr = ( imag(a.u11) + imag(a.u22) + 2.0*imag(a.u33) )/3.0 ditr = ( imag(a.u11) + imag(a.u22) - 2.0*imag(a.u33) )/3.0
m12 = (a.u12 - conj(a.u21))/2.0 m12 = (a.u12 - conj(a.u21))/2.0
m13 = (a.u13 - conj(a.u31))/2.0 m13 = (a.u13 - conj(a.u31))/2.0
m23 = (a.u23 - conj(a.u32))/2.0 m23 = (a.u23 - conj(a.u32))/2.0

1
src/LatticeGPU.jl
View file

@ -38,6 +38,7 @@ export omf4, leapfrog, omf2
include("YM/YM.jl") include("YM/YM.jl")
using .YM using .YM
export ztwist
export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2 export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2
export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3, Eoft_clover, Eoft_plaq, Qtop export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3, Eoft_clover, Eoft_plaq, Qtop

41
src/Space/Space.jl
View file

@ -33,11 +33,17 @@ struct SpaceParm{N,M,B,D}
bsz::Int64 bsz::Int64
rsz::Int64 rsz::Int64
function SpaceParm{N}(x, y) where {N} ntw::NTuple{M,Int64}
function SpaceParm{N}(x, y, nt::Union{Nothing,NTuple{I,Int64}}=nothing) where {N,I}
M = convert(Int64, round(N*(N-1)/2)) M = convert(Int64, round(N*(N-1)/2))
N == length(x) || throw(ArgumentError("Lattice size incorrect length for dimension $N")) N == length(x) || throw(ArgumentError("Lattice size incorrect length for dimension $N"))
N == length(y) || throw(ArgumentError("Block size incorrect length for dimension $N")) N == length(y) || throw(ArgumentError("Block size incorrect length for dimension $N"))
if any(i->i!=0, x.%y)
error("Lattice size not divisible by block size.")
end
pls = Vector{Tuple{Int64, Int64}}() pls = Vector{Tuple{Int64, Int64}}()
for i in N:-1:1 for i in N:-1:1
for j in 1:i-1 for j in 1:i-1
@ -56,15 +62,32 @@ struct SpaceParm{N,M,B,D}
end end
D = prod(y) D = prod(y)
if nt == nothing
ntw = ntuple(i->0, M)
else
ntw = nt
end
return new{N,M,BC_PERIODIC,D}(N, x, M, tuple(pls...), y, return new{N,M,BC_PERIODIC,D}(N, x, M, tuple(pls...), y,
tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r)) tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r), ntw)
end end
function SpaceParm{N}(x, y, ibc) where {N} function SpaceParm{N}(x, y, ibc::Int, nt::Union{Nothing,NTuple{I,Int64}}=nothing) where {N,I}
M = convert(Int64, round(N*(N-1)/2)) M = convert(Int64, round(N*(N-1)/2))
N == length(x) || throw(ArgumentError("Lattice size incorrect length for dimension $N")) N == length(x) || throw(ArgumentError("Lattice size incorrect length for dimension $N"))
N == length(y) || throw(ArgumentError("Block size incorrect length for dimension $N")) N == length(y) || throw(ArgumentError("Block size incorrect length for dimension $N"))
if any(i->i!=0, x.%y)
error("Lattice size not divisible by block size.")
end
if nt!=nothing
if (ibc==BC_SF_AFWB) || (ibc==BC_SF_ORBI)
if any(i->i!=0, nt[1:N-1])
error("Planes in T direction cannot be twisted with SF boundary conditions")
end
end
end
pls = Vector{Tuple{Int64, Int64}}() pls = Vector{Tuple{Int64, Int64}}()
for i in N:-1:1 for i in N:-1:1
for j in 1:i-1 for j in 1:i-1
@ -83,11 +106,17 @@ struct SpaceParm{N,M,B,D}
end end
D = prod(y) D = prod(y)
if nt == nothing
ntw = ntuple(i->0,M)
else
ntw = nt
end
return new{N,M,ibc,D}(N, x, M, tuple(pls...), y, return new{N,M,ibc,D}(N, x, M, tuple(pls...), y,
tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r)) tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r), ntw)
end end
end end
export SpaceParm export SpaceParm
function Base.show(io::IO, lp::SpaceParm{N,M,B,D}) where {N,M,B,D} function Base.show(io::IO, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}
println(io, "Lattice dimensions: ", lp.ndim) println(io, "Lattice dimensions: ", lp.ndim)
@ -113,6 +142,8 @@ function Base.show(io::IO, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}
end end
println(io,lp.blk[end], " [", lp.bsz, println(io,lp.blk[end], " [", lp.bsz,
"] (Number of blocks: [", lp.rsz,"])") "] (Number of blocks: [", lp.rsz,"])")
println(io, "Twist tensor: ", lp.ntw)
return return
end end
@ -170,7 +201,7 @@ Given a point `x` with index `p`, this routine returns the index of the point
end end
""" """
function up(p::NTuple{2,Int64}, id::Int64, lp::SpaceParm) function updw(p::NTuple{2,Int64}, id::Int64, lp::SpaceParm)
Given a point `x` with index `p`, this routine returns the index of the points Given a point `x` with index `p`, this routine returns the index of the points
`x + a id` and `x - a id`. `x + a id` and `x - a id`.

24
src/YM/YM.jl
View file

@ -104,15 +104,29 @@ function Base.show(io::IO, ymws::YMworkspace)
println(io, "Workspace for Group: ", ymws.GRP) println(io, "Workspace for Group: ", ymws.GRP)
println(io, " Algebra: ", ymws.ALG) println(io, " Algebra: ", ymws.ALG)
println(io, "Precision: ", ymws.PRC) println(io, "Precision: ", ymws.PRC)
if ymws.fpln == nothing
println(io, " - Running in memory efficient mode")
else
println(io, " - Running in computing efficient mode")
end
return nothing return nothing
end end
function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}) where {T,G,N,M,B,D}
function plnf(ipl)
id1, id2 = lp.plidx[ipl]
return convert(Complex{T},exp(2im * pi * lp.ntw[ipl]/(lp.iL[id1]*lp.iL[id2]*gp.ng)))
end
return ntuple(i->plnf(i), M)
end
function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}, ipl::Int) where {T,G,N,M,B,D}
id1, id2 = lp.plidx[ipl]
return convert(Complex{T},exp(2im * pi * lp.ntw[ipl]/(lp.iL[id1]*lp.iL[id2]*gp.ng)))
end
export ztwist
include("YMfields.jl") include("YMfields.jl")
export randomize!, zero!, norm2 export randomize!, zero!, norm2

65
src/YM/YMact.jl
View file

@ -9,13 +9,14 @@
### created: Mon Jul 12 18:31:19 2021 ### created: Mon Jul 12 18:31:19 2021
### ###
function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ztw, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
it = point_time((b, r), lp) it = point_time((b, r), lp)
Ush = @cuStaticSharedMem(T, (D,2)) Ush = @cuStaticSharedMem(T, (D,2))
ipl = 0
S = zero(eltype(plx)) S = zero(eltype(plx))
for id1 in N:-1:1 for id1 in N:-1:1
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
@ -26,6 +27,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
Ush[b,2] = U[b,id2,r] Ush[b,2] = U[b,id2,r]
sync_threads() sync_threads()
ipl = ipl + 1
# H2 staple # H2 staple
(b1, r1) = up((b,r), id1, lp) (b1, r1) = up((b,r), id1, lp)
@ -100,7 +102,8 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm
elseif (it == 1) && SFBC elseif (it == 1) && SFBC
S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3)) + c1*tr(g2*h2/gb) S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3)) + c1*tr(g2*h2/gb)
else else
S += c0*tr(g2*ga/gb) + c1*( tr(g2*h2/gb) + tr(g2*ga/h3)) S += ztw[ipl]*c0*tr(g2*ga/gb) +
(ztw[ipl]^2*c1)*( tr(g2*h2/gb) + tr(g2*ga/h3))
end end
end end
@ -112,14 +115,14 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm
return nothing return nothing
end end
function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, ztw, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
it = point_time((b, r), lp)
Ush = @cuStaticSharedMem(T, (D,2)) Ush = @cuStaticSharedMem(T, (D,2))
S = zero(eltype(plx)) S = zero(eltype(plx))
ipl = 0
for id1 in N:-1:1 for id1 in N:-1:1
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
Ush[b,1] = U[b,id1,r] Ush[b,1] = U[b,id1,r]
@ -131,6 +134,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
Ush[b,2] = U[b,id2,r] Ush[b,2] = U[b,id2,r]
sync_threads() sync_threads()
ipl = ipl + 1
if ru1 == r if ru1 == r
gt1 = Ush[bu1,2] gt1 = Ush[bu1,2]
@ -150,7 +154,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
if SFBND if SFBND
S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
else else
S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) S += ztw[ipl]*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
end end
end end
end end
@ -161,7 +165,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
return nothing return nothing
end end
function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
it = point_time((b,r), lp) it = point_time((b,r), lp)
@ -198,32 +202,31 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
g2 = Ush[b,2]\Ush[b,1] g2 = Ush[b,2]\Ush[b,1]
if SFBC && (it == 1) if SFBC && (it == 1)
X = cG*projalg(Ush[b,1]*g1/Ush[b,2]) X = cG*projalg(ztw,Ush[b,1]*g1/Ush[b,2])
frc1[b ,id1, r ] -= X frc1[b ,id1, r ] -= X
frc2[bu1,id2,ru1] -= cG*projalg(g1*g2) frc2[bu1,id2,ru1] -= cG*projalg(ztw,g1*g2)
frc2[bu2,id1,ru2] += cG*projalg(g2*g1) frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1)
elseif SFBC && (it == lp.iL[end]) elseif SFBC && (it == lp.iL[end])
X = cG*projalg(Ush[b,1]*g1/Ush[b,2]) X = cG*projalg(ztw,Ush[b,1]*g1/Ush[b,2])
frc1[b ,id1, r ] -= X frc1[b ,id1, r ] -= X
frc1[b ,id2, r ] += X frc1[b ,id2, r ] += X
frc2[bu2,id1,ru2] += cG*projalg(g2*g1) frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1)
else else
X = projalg(Ush[b,1]*g1/Ush[b,2]) X = projalg(ztw,Ush[b,1]*g1/Ush[b,2])
frc1[b ,id1, r ] -= X frc1[b ,id1, r ] -= X
frc1[b ,id2, r ] += X frc1[b ,id2, r ] += X
frc2[bu1,id2,ru1] -= projalg(g1*g2) frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2)
frc2[bu2,id1,ru2] += projalg(g2*g1) frc2[bu2,id1,ru2] += projalg(ztw,g2*g1)
end end
end end
return nothing return nothing
end end
function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
it = point_time((b, r), lp) it = point_time((b, r), lp)
@ -371,19 +374,18 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T,
frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((Ush[b,2]\h1)*g1) + frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((Ush[b,2]\h1)*g1) +
c1 * projalg(h4\Ush[b,1]*g1) c1 * projalg(h4\Ush[b,1]*g1)
else else
X = c0*projalg(Ush[b,1]*g1/Ush[b,2]) + c1 * (projalg(Ush[b,1]*h2/(Ush[b,2]*gb)) + zsq = ztw[ipl]^2
projalg(Ush[b,1]*ga/(Ush[b,2]*h3))) X = projalg(c0*ztw[ipl],Ush[b,1]*g1/Ush[b,2]) + projalg(zsq*c1,Ush[b,1]*h2/(Ush[b,2]*gb)) +
projalg(zsq*c1,Ush[b,1]*ga/(Ush[b,2]*h3))
frc1[b,id1,r] -= X + c1*projalg(Ush[b,1]*g1/h4) frc1[b,id1,r] -= X + projalg(zsq*c1,Ush[b,1]*g1/h4)
frc1[b,id2,r] += X + c1*projalg(h1*g1/Ush[b,2]) frc1[b,id2,r] += X + projalg(zsq*c1,h1*g1/Ush[b,2])
frc2[bu1,id2,ru1] -= c0*projalg(g1*g2) + c1*( projalg((ga/h3)*g2) + frc2[bu1,id2,ru1] -= projalg(c0*ztw[ipl],g1*g2) + projalg(zsq*c1,(ga/h3)*g2) +
projalg((g1/h4)*Ush[b,1]) + projalg(zsq*c1,(g1/h4)*Ush[b,1]) + projalg(zsq*c1,(g1/Ush[b,2])*h1)
projalg((g1/Ush[b,2])*h1) )
frc2[bu2,id1,ru2] += c0*projalg(g2*g1) + c1*( projalg((Ush[b,2]\h1)*g1) + frc2[bu2,id1,ru2] += projalg(c0*ztw[ipl],g2*g1) + projalg(zsq*c1,(Ush[b,2]\h1)*g1) +
projalg(g2*h2/gb) + projalg(zsq*c1,g2*h2/gb) + projalg(zsq*c1,h4\Ush[b,1]*g1)
projalg(h4\Ush[b,1]*g1) )
end end
end end
@ -399,13 +401,14 @@ the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1`
""" """
function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm) function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm)
ztw = ztwist(gp, lp)
if abs(c0-1) < 1.0E-10 if abs(c0-1) < 1.0E-10
@timeit "Wilson gauge force" begin @timeit "Wilson gauge force" begin
force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, lp::SpaceParm) force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm)
end end
else else
@timeit "Improved gauge force" begin @timeit "Improved gauge force" begin
force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, lp::SpaceParm, c0) force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm, c0)
end end
end end
return nothing return nothing
@ -416,20 +419,20 @@ force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm) = force_gauge(y
force_wilson(ymws::YMworkspace, U, cG, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp.cG[1], gp, lp) force_wilson(ymws::YMworkspace, U, cG, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp.cG[1], gp, lp)
function force_pln!(frc1, ftmp, U, Ubnd, cG, lp::SpaceParm, c0=1) function force_pln!(frc1, ftmp, U, Ubnd, cG, ztw, lp::SpaceParm, c0=1)
fill!(frc1, zero(eltype(frc1))) fill!(frc1, zero(eltype(frc1)))
fill!(ftmp, zero(eltype(ftmp))) fill!(ftmp, zero(eltype(ftmp)))
if c0 == 1 if c0 == 1
for i in 1:lp.npls for i in 1:lp.npls
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_pln!(frc1,ftmp,U, Ubnd, cG,i,lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_pln!(frc1,ftmp,U, Ubnd, cG, ztw[i], i,lp)
end end
end end
else else
for i in 1:lp.npls for i in 1:lp.npls
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_impr_pln!(frc1,ftmp,U,c0,(1-c0)/8,Ubnd, cG, i,lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_impr_pln!(frc1,ftmp,U,c0,(1-c0)/8,Ubnd, cG, ztw[i], i,lp)
end end
end end
end end

89
src/YM/YMflow.jl
View file

@ -146,6 +146,9 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMwo
@timeit "E(t) plaquette measurement" begin @timeit "E(t) plaquette measurement" begin
ztw = ztwist(gp, lp)
SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) )
tp = ntuple(i->i, N-1) tp = ntuple(i->i, N-1)
V3 = prod(lp.iL[1:end-1]) V3 = prod(lp.iL[1:end-1])
@ -154,7 +157,7 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMwo
for ipl in 1:M for ipl in 1:M
fill!(Etmp, zero(T)) fill!(Etmp, zero(T))
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq_pln!(ymws.cm, U, ipl, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq_pln!(ymws.cm, U, gp.Ubnd, ztw[ipl], ipl, lp)
end end
Etmp .= (gp.ng .- reshape(Array(CUDA.mapreduce(real, +, ymws.cm;dims=tp)),lp.iL[end])/V3 ) Etmp .= (gp.ng .- reshape(Array(CUDA.mapreduce(real, +, ymws.cm;dims=tp)),lp.iL[end])/V3 )
@ -162,7 +165,9 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMwo
for it in 2:lp.iL[end] for it in 2:lp.iL[end]
Eslc[it,ipl] = Etmp[it] + Etmp[it-1] Eslc[it,ipl] = Etmp[it] + Etmp[it-1]
end end
if !SFBC
Eslc[1,ipl] = Etmp[1] + Etmp[end] Eslc[1,ipl] = Etmp[1] + Etmp[end]
end
else else
for it in 1:lp.iL[end] for it in 1:lp.iL[end]
Eslc[it,ipl] = 2*Etmp[it] Eslc[it,ipl] = 2*Etmp[it]
@ -179,17 +184,24 @@ end
Eoft_plaq(U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) where {T,N,M,B,D} = Eoft_plaq(zeros(T,lp.iL[end],M), U, gp, lp, ymws) Eoft_plaq(U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) where {T,N,M,B,D} = Eoft_plaq(zeros(T,lp.iL[end],M), U, gp, lp, ymws)
function krnl_plaq_pln!(plx, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd::T, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
id1, id2 = lp.plidx[ipl] id1, id2 = lp.plidx[ipl]
SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI)) && (id1 == lp.iL[end])
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
if SFBC && (ru1 != r)
gt = Ubnd
else
gt = U[bu1,id2,ru1]
end
I = point_coord((b,r), lp) I = point_coord((b,r), lp)
plx[I] = tr(U[b,1,r]*U[bu1,id2,ru1] / (U[b,2,r]*U[bu2,id1,ru2])) plx[I] = ztw*tr(U[b,1,r]*gt / (U[b,2,r]*U[bu2,id1,ru2]))
return nothing return nothing
end end
@ -200,39 +212,40 @@ end
Measure the topological charge `Q` of the configuration `U`. If the argument `Qslc` is present Measure the topological charge `Q` of the configuration `U`. If the argument `Qslc` is present
the contribution for each Euclidean time slice are returned. the contribution for each Euclidean time slice are returned.
""" """
function Qtop(Qslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace) where {M,B,D} function Qtop(Qslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace) where {M,B,D}
@timeit "Qtop measurement" begin @timeit "Qtop measurement" begin
ztw = ztwist(gp, lp)
tp = (1,2,3) tp = (1,2,3)
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 1,5, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,6, ztw[1], ztw[5], lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, +, ymws.frc1, ymws.frc2, U, lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 2,4, lp)
end end
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp)
end end
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 3,6, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 2,5, ztw[2], ztw[4], lp)
end end
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, +, ymws.frc1, ymws.frc2, U, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, +, ymws.frc1, ymws.frc2, U, lp)
end end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,4, ztw[3], ztw[6], lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp)
end
Qslc .= reshape(Array(CUDA.reduce(+, ymws.rm; dims=tp)),lp.iL[end])./(32*pi^2) Qslc .= reshape(Array(CUDA.reduce(+, ymws.rm; dims=tp)),lp.iL[end])./(32*pi^2)
end end
return sum(Qslc) return sum(Qslc)
end end
Qtop(U, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where {T,M,D} = Qtop(zeros(T,lp.iL[end],M), U, lp, ymws) Qtop(U, gp::GaugeParm, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where {T,M,D} = Qtop(zeros(T,lp.iL[end],M), U, gp, lp, ymws)
""" """
@ -241,7 +254,7 @@ Qtop(U, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where {T,M,D} = Qtop(zeros(T
Measure the action density `E(t)` using the clover discretization. If the argument `Eslc` Measure the action density `E(t)` using the clover discretization. If the argument `Eslc`
the contribution for each Euclidean time slice and plane are returned. the contribution for each Euclidean time slice and plane are returned.
""" """
function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace{T}) where {T,M,B,D} function Eoft_clover(Eslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace{T}) where {T,M,B,D}
function acum(ipl1, ipl2, Etmp) function acum(ipl1, ipl2, Etmp)
@ -270,22 +283,23 @@ function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace{T}) wher
@timeit "E(t) clover measurement" begin @timeit "E(t) clover measurement" begin
ztw = ztwist(gp, lp)
fill!(Eslc,zero(T)) fill!(Eslc,zero(T))
Etmp = zeros(T,lp.iL[end]) Etmp = zeros(T,lp.iL[end])
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 1,2, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,2, ztw[1], ztw[2], lp)
end end
acum(1,2,Etmp) acum(1,2,Etmp)
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 3,4, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,4, ztw[3], ztw[4], lp)
end end
acum(3,4,Etmp) acum(3,4,Etmp)
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 5,6, lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 5,6, ztw[5], ztw[6], lp)
end end
acum(5,6,Etmp) acum(5,6,Etmp)
@ -293,7 +307,7 @@ function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace{T}) wher
return sum(Eslc)/lp.iL[end] return sum(Eslc)/lp.iL[end]
end end
Eoft_clover(U, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace{T}) where {T,N,M,B,D} = Eoft_clover(zeros(T,lp.iL[end],M), U, lp, ymws) Eoft_clover(U, gp::GaugeParm, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace{T}) where {T,N,M,B,D} = Eoft_clover(zeros(T,lp.iL[end],M), U, gp, lp, ymws)
function krnl_add_et!(rm, op, frc1, U, lp::SpaceParm{4,M,B,D}) where {M,B,D} function krnl_add_et!(rm, op, frc1, U, lp::SpaceParm{4,M,B,D}) where {M,B,D}
@ -312,14 +326,16 @@ function krnl_add_qd!(rm, op, frc1, frc2, U, lp::SpaceParm{4,M,B,D}) where {M,B,
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
I = point_coord((b,r), lp) I = point_coord((b,r), lp)
rm[I] = op(dot( (frc1[b,1,r]+frc1[b,2,r]+frc1[b,3,r]+frc1[b,4,r]), rm[I] += op(dot( (frc1[b,1,r]+frc1[b,2,r]+frc1[b,3,r]+frc1[b,4,r]),
(frc2[b,1,r]+frc2[b,2,r]+frc2[b,3,r]+frc2[b,4,r]) ) ) (frc2[b,1,r]+frc2[b,2,r]+frc2[b,3,r]+frc2[b,4,r]) ) )
return nothing return nothing
end end
function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::SpaceParm{4,M,B,D}) where {T,M,B,D} function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, ipl1, ipl2, ztw1, ztw2, lp::SpaceParm{4,M,B,D}) where {T,M,B,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
it = point_time((b,r), lp)
SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) )
Ush = @cuStaticSharedMem(T, (D,2)) Ush = @cuStaticSharedMem(T, (D,2))
@ -329,27 +345,39 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::Spa
Ush[b,2] = U[b,id2,r] Ush[b,2] = U[b,id2,r]
sync_threads() sync_threads()
SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N)
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
bd, rd = up((bu1, ru1), id2, lp) bd, rd = up((bu1, ru1), id2, lp)
if ru1 == r if ru1 == r
gt1 = Ush[bu1,2] gt1 = Ush[bu1,2]
else else
gt1 = U[bu1,id2,ru1] gt1 = U[bu1,id2,ru1]
end end
if ru2 == r if ru2 == r
gt2 = Ush[bu2,1] gt2 = Ush[bu2,1]
else
if SFBC && (it == lp.iL[end])
gt2 = Ubnd
else else
gt2 = U[bu2,id1,ru2] gt2 = U[bu2,id1,ru2]
end end
end
l1 = gt1/gt2 l1 = gt1/gt2
l2 = Ush[b,2]\Ush[b,1] l2 = Ush[b,2]\Ush[b,1]
if SFBC && (it == lp.iL[end])
frc1[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) frc1[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2])
frc1[bu1,2,ru1] = projalg(l1*l2)
frc1[bd,3,rd] = projalg(gt2\(l2*gt1))
frc1[bu2,4,ru2] = projalg(l2*l1) frc1[bu2,4,ru2] = projalg(l2*l1)
else
frc1[b,1,r] = projalg(ztw1, Ush[b,1]*l1/Ush[b,2])
frc1[bu1,2,ru1] = projalg(ztw1, l1*l2)
frc1[bd,3,rd] = projalg(ztw1, gt2\(l2*gt1))
frc1[bu2,4,ru2] = projalg(ztw1, l2*l1)
end
# Second plane # Second plane
id1, id2 = lp.plidx[ipl2] id1, id2 = lp.plidx[ipl2]
@ -357,6 +385,8 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::Spa
Ush[b,2] = U[b,id2,r] Ush[b,2] = U[b,id2,r]
sync_threads() sync_threads()
SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N)
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
bd, rd = up((bu1, ru1), id2, lp) bd, rd = up((bu1, ru1), id2, lp)
@ -367,17 +397,26 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::Spa
end end
if ru2 == r if ru2 == r
gt2 = Ush[bu2,1] gt2 = Ush[bu2,1]
else
if SFBC && (it == lp.iL[end])
gt2 = Ubnd
else else
gt2 = U[bu2,id1,ru2] gt2 = U[bu2,id1,ru2]
end end
end
l1 = gt1/gt2 l1 = gt1/gt2
l2 = Ush[b,2]\Ush[b,1] l2 = Ush[b,2]\Ush[b,1]
if SFBC && (it == lp.iL[end])
frc2[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) frc2[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2])
frc2[bu1,2,ru1] = projalg(l1*l2)
frc2[bd,3,rd] = projalg(gt2\(l2*gt1))
frc2[bu2,4,ru2] = projalg(l2*l1) frc2[bu2,4,ru2] = projalg(l2*l1)
else
frc2[b,1,r] = projalg(ztw1, Ush[b,1]*l1/Ush[b,2])
frc2[bu1,2,ru1] = projalg(ztw1, l1*l2)
frc2[bd,3,rd] = projalg(ztw1, gt2\(l2*gt1))
frc2[bu2,4,ru2] = projalg(ztw1, l2*l1)
end
return nothing return nothing
end end

8
src/YM/YMhmc.jl
View file

@ -17,16 +17,17 @@ Returns the value of the gauge plaquette action for the configuration U. The par
""" """
function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat
ztw = ztwist(gp, lp)
if abs(gp.c0-1) < 1.0E-10 if abs(gp.c0-1) < 1.0E-10
@timeit "Wilson gauge action" begin @timeit "Wilson gauge action" begin
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, gp.cG[1], lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, gp.cG[1], ztw, lp)
end end
end end
else else
@timeit "Improved gauge action" begin @timeit "Improved gauge action" begin
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_impr!(ymws.cm, U, gp.c0, (1-gp.c0)/8, gp.Ubnd, gp.cG[1], lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_impr!(ymws.cm, U, gp.c0, (1-gp.c0)/8, gp.Ubnd, gp.cG[1], ztw, lp)
end end
end end
end end
@ -38,9 +39,10 @@ end
function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace)
ztw = ztwist(gp, lp)
@timeit "Plaquette measurement" begin @timeit "Plaquette measurement" begin
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, one(gp.cG[1]), lp) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, one(gp.cG[1]), ztw, lp)
end end
end end

32
src/main/times.jl
View file

@ -7,7 +7,8 @@ Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
using LatticeGPU using LatticeGPU
# Set lattice/block size # Set lattice/block size
lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_SF_AFWB) ntwist = (0,0,0,1,0,0)
lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist)
println("Space Parameters: ", lp) println("Space Parameters: ", lp)
# Seed RNG # Seed RNG
@ -21,6 +22,10 @@ ALG = SU3alg
PREC = Float64 PREC = Float64
println("Precision: ", PREC) println("Precision: ", PREC)
# MD integrator
int = omf4(PREC, 0.1, 10)
println(int)
println("Allocating YM workspace") println("Allocating YM workspace")
ymws = YMworkspace(GRP, PREC, lp) ymws = YMworkspace(GRP, PREC, lp)
@ -44,19 +49,17 @@ println("Initial Action: ")
@time S = gauge_action(U, lp, gp, ymws) @time S = gauge_action(U, lp, gp, ymws)
dt = 0.1 HMC!(U,int.eps,1,lp, gp, ymws, noacc=true)
ns = 10
HMC!(U,dt,1,lp, gp, ymws, noacc=true)
pl = Vector{Float64}() pl = Vector{Float64}()
for i in 1:4 for i in 1:4
@time dh, acc = HMC!(U,dt,ns,lp, gp, ymws, noacc=true) @time dh, acc = HMC!(U,int,lp, gp, ymws, noacc=true)
println("# HMC: ", acc, " ", dh) println("# HMC: ", acc, " ", dh)
push!(pl, plaquette(U,lp, gp, ymws)) push!(pl, plaquette(U,lp, gp, ymws))
println("# Plaquette: ", pl[end], "\n") println("# Plaquette: ", pl[end], "\n")
@time x, y = sfcoupling(U,lp,gp,ymws) # @time x, y = sfcoupling(U,lp,gp,ymws)
println("SF coupling: ", x, " ", y) # println("SF coupling: ", x, " ", y)
end end
wfl_rk3(U, 1, 0.01, gp, lp, ymws) wfl_rk3(U, 1, 0.01, gp, lp, ymws)
@ -65,14 +68,14 @@ println("Action: ", gauge_action(U, lp, gp, ymws))
println("Time for 100 steps of RK3 flow integrator: ") println("Time for 100 steps of RK3 flow integrator: ")
@time wfl_rk3(U, 100, 0.01, gp, lp, ymws) @time wfl_rk3(U, 100, 0.01, gp, lp, ymws)
eoft = Eoft_plaq(U, gp, lp, ymws) eoft = Eoft_plaq(U, gp, lp, ymws)
eoft = Eoft_clover(U, lp, ymws) eoft = Eoft_clover(U, gp, lp, ymws)
qtop = Qtop(U, lp, ymws) qtop = Qtop(U, gp, lp, ymws)
@time eoft = Eoft_plaq(U, gp, lp, ymws) @time eoft = Eoft_plaq(U, gp, lp, ymws)
println("Plaq: ", eoft) println("Plaq: ", eoft)
@time eoft = Eoft_clover(U, lp, ymws) @time eoft = Eoft_clover(U, gp, lp, ymws)
println("Clov: ", eoft) println("Clov: ", eoft)
@time qtop = Qtop(U, lp, ymws) @time qtop = Qtop(U, gp, lp, ymws)
println("Qtop: ", qtop) println("Qtop: ", qtop)
println("Action: ", gauge_action(U, lp, gp, ymws)) println("Action: ", gauge_action(U, lp, gp, ymws))
@ -88,13 +91,10 @@ println("Initial Action: ")
@time S = gauge_action(U, lp, gp, ymws) @time S = gauge_action(U, lp, gp, ymws)
dt = 0.1 HMC!(U,int.eps,1,lp, gp, ymws, noacc=true)
ns = 10
HMC!(U,dt,1,lp, gp, ymws, noacc=true)
pl = Vector{Float64}() pl = Vector{Float64}()
for i in 1:4 for i in 1:4
@time dh, acc = HMC!(U,dt,ns,lp, gp, ymws, noacc=true) @time dh, acc = HMC!(U,int,lp, gp, ymws, noacc=true)
println("# HMC: ", acc, " ", dh) println("# HMC: ", acc, " ", dh)
push!(pl, plaquette(U,lp, gp, ymws)) push!(pl, plaquette(U,lp, gp, ymws))
println("# Plaquette: ", pl[end], "\n") println("# Plaquette: ", pl[end], "\n")