From 410a6d7b25f273de49988dcb85dbca9a506538f1 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Sun, 31 Oct 2021 08:37:53 +0100 Subject: [PATCH] Full support for twisted bc --- src/Groups/AlgebraSU3.jl | 24 +++++++++ src/Groups/M3x3.jl | 2 +- src/LatticeGPU.jl | 1 + src/Space/Space.jl | 41 +++++++++++++-- src/YM/YM.jl | 24 +++++++-- src/YM/YMact.jl | 73 +++++++++++++------------- src/YM/YMflow.jl | 109 ++++++++++++++++++++++++++------------- src/YM/YMhmc.jl | 8 +-- src/main/times.jl | 32 ++++++------ 9 files changed, 214 insertions(+), 100 deletions(-) diff --git a/src/Groups/AlgebraSU3.jl b/src/Groups/AlgebraSU3.jl index 107ab4a..52e8005 100644 --- a/src/Groups/AlgebraSU3.jl +++ b/src/Groups/AlgebraSU3.jl @@ -24,6 +24,30 @@ function projalg(a::SU3{T}) where T <: AbstractFloat sr3ov2*(ditr)) 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 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) diff --git a/src/Groups/M3x3.jl b/src/Groups/M3x3.jl index 30866f0..32850bc 100644 --- a/src/Groups/M3x3.jl +++ b/src/Groups/M3x3.jl @@ -117,7 +117,7 @@ function projalg(a::M3x3{T}) where T <: AbstractFloat 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 m13 = (a.u13 - conj(a.u31))/2.0 m23 = (a.u23 - conj(a.u32))/2.0 diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 9b9f31a..a86f9f9 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -38,6 +38,7 @@ export omf4, leapfrog, omf2 include("YM/YM.jl") using .YM +export ztwist export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2 export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3, Eoft_clover, Eoft_plaq, Qtop diff --git a/src/Space/Space.jl b/src/Space/Space.jl index f35a603..5ff9001 100644 --- a/src/Space/Space.jl +++ b/src/Space/Space.jl @@ -32,12 +32,18 @@ struct SpaceParm{N,M,B,D} bsz::Int64 rsz::Int64 + + ntw::NTuple{M,Int64} - function SpaceParm{N}(x, y) where {N} + function SpaceParm{N}(x, y, nt::Union{Nothing,NTuple{I,Int64}}=nothing) where {N,I} M = convert(Int64, round(N*(N-1)/2)) N == length(x) || throw(ArgumentError("Lattice 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}}() for i in N:-1:1 for j in 1:i-1 @@ -56,15 +62,32 @@ struct SpaceParm{N,M,B,D} end 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, - tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r)) + tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r), ntw) 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)) N == length(x) || throw(ArgumentError("Lattice 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}}() for i in N:-1:1 for j in 1:i-1 @@ -83,11 +106,17 @@ struct SpaceParm{N,M,B,D} end 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, - tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r)) + tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r), ntw) end end export SpaceParm + function Base.show(io::IO, lp::SpaceParm{N,M,B,D}) where {N,M,B,D} 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 println(io,lp.blk[end], " [", lp.bsz, "] (Number of blocks: [", lp.rsz,"])") + println(io, "Twist tensor: ", lp.ntw) + return end @@ -170,7 +201,7 @@ Given a point `x` with index `p`, this routine returns the index of the point 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 `x + a id` and `x - a id`. diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 27c4d58..4f11846 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -104,15 +104,29 @@ function Base.show(io::IO, ymws::YMworkspace) println(io, "Workspace for Group: ", ymws.GRP) println(io, " Algebra: ", ymws.ALG) 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 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") export randomize!, zero!, norm2 diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 8672791..a23fb5a 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -9,13 +9,14 @@ ### 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 it = point_time((b, r), lp) Ush = @cuStaticSharedMem(T, (D,2)) - + + ipl = 0 S = zero(eltype(plx)) for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) @@ -26,7 +27,8 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm bu2, ru2 = up((b, r), id2, lp) Ush[b,2] = U[b,id2,r] sync_threads() - + ipl = ipl + 1 + # H2 staple (b1, r1) = up((b,r), id1, lp) if r1 == r @@ -100,7 +102,8 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm elseif (it == 1) && SFBC S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3)) + c1*tr(g2*h2/gb) 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 @@ -112,14 +115,14 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, lp::SpaceParm return nothing 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 - it = point_time((b, r), lp) Ush = @cuStaticSharedMem(T, (D,2)) S = zero(eltype(plx)) + ipl = 0 for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) Ush[b,1] = U[b,id1,r] @@ -131,7 +134,8 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D bu2, ru2 = up((b, r), id2, lp) Ush[b,2] = U[b,id2,r] sync_threads() - + ipl = ipl + 1 + if ru1 == r gt1 = Ush[bu1,2] else @@ -150,7 +154,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D if SFBND S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) 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 @@ -161,10 +165,10 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D return nothing 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 - it = point_time((b, r), lp) + it = point_time((b,r), lp) Ush = @cuStaticSharedMem(T, (D,2)) @@ -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] 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 - frc2[bu1,id2,ru1] -= cG*projalg(g1*g2) - frc2[bu2,id1,ru2] += cG*projalg(g2*g1) + frc2[bu1,id2,ru1] -= cG*projalg(ztw,g1*g2) + frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) 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 ,id2, r ] += X - frc2[bu2,id1,ru2] += cG*projalg(g2*g1) + frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) 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 ,id2, r ] += X - frc2[bu1,id2,ru1] -= projalg(g1*g2) - frc2[bu2,id1,ru2] += projalg(g2*g1) + frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2) + frc2[bu2,id1,ru2] += projalg(ztw,g2*g1) end - end return nothing 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 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) + c1 * projalg(h4\Ush[b,1]*g1) else - X = c0*projalg(Ush[b,1]*g1/Ush[b,2]) + c1 * (projalg(Ush[b,1]*h2/(Ush[b,2]*gb)) + - projalg(Ush[b,1]*ga/(Ush[b,2]*h3))) + zsq = ztw[ipl]^2 + 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,id2,r] += X + c1*projalg(h1*g1/Ush[b,2]) + frc1[b,id1,r] -= X + projalg(zsq*c1,Ush[b,1]*g1/h4) + 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) + - projalg((g1/h4)*Ush[b,1]) + - projalg((g1/Ush[b,2])*h1) ) + frc2[bu1,id2,ru1] -= projalg(c0*ztw[ipl],g1*g2) + projalg(zsq*c1,(ga/h3)*g2) + + projalg(zsq*c1,(g1/h4)*Ush[b,1]) + projalg(zsq*c1,(g1/Ush[b,2])*h1) - frc2[bu2,id1,ru2] += c0*projalg(g2*g1) + c1*( projalg((Ush[b,2]\h1)*g1) + - projalg(g2*h2/gb) + - projalg(h4\Ush[b,1]*g1) ) + frc2[bu2,id1,ru2] += projalg(c0*ztw[ipl],g2*g1) + projalg(zsq*c1,(Ush[b,2]\h1)*g1) + + projalg(zsq*c1,g2*h2/gb) + projalg(zsq*c1,h4\Ush[b,1]*g1) 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) + ztw = ztwist(gp, lp) if abs(c0-1) < 1.0E-10 @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 else @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 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) -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!(ftmp, zero(eltype(ftmp))) if c0 == 1 for i in 1:lp.npls 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 else for i in 1:lp.npls 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 diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 9a6a24e..c777a3f 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -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 + ztw = ztwist(gp, lp) + SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) + tp = ntuple(i->i, N-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 fill!(Etmp, zero(T)) 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 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] Eslc[it,ipl] = Etmp[it] + Etmp[it-1] end - Eslc[1,ipl] = Etmp[1] + Etmp[end] + if !SFBC + Eslc[1,ipl] = Etmp[1] + Etmp[end] + end else for it in 1:lp.iL[end] 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) -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 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) 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) - 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 end @@ -200,39 +212,40 @@ end Measure the topological charge `Q` of the configuration `U`. If the argument `Qslc` is present 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 + ztw = ztwist(gp, lp) tp = (1,2,3) CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, 1,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) + 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, 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 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, 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) end return sum(Qslc) 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` 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) @@ -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 + ztw = ztwist(gp, lp) fill!(Eslc,zero(T)) Etmp = zeros(T,lp.iL[end]) 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 acum(1,2,Etmp) 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 acum(3,4,Etmp) 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 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] 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} @@ -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 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]), - (frc2[b,1,r]+frc2[b,2,r]+frc2[b,3,r]+frc2[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]) ) ) return nothing 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 + it = point_time((b,r), lp) + SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) Ush = @cuStaticSharedMem(T, (D,2)) @@ -329,34 +345,48 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::Spa Ush[b,2] = U[b,id2,r] sync_threads() + SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N) + bu1, ru1 = up((b, r), id1, lp) bu2, ru2 = up((b, r), id2, lp) bd, rd = up((bu1, ru1), id2, lp) if ru1 == r gt1 = Ush[bu1,2] else + gt1 = U[bu1,id2,ru1] end if ru2 == r gt2 = Ush[bu2,1] else - gt2 = U[bu2,id1,ru2] + if SFBC && (it == lp.iL[end]) + gt2 = Ubnd + else + gt2 = U[bu2,id1,ru2] + end end l1 = gt1/gt2 l2 = Ush[b,2]\Ush[b,1] - 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) - + if SFBC && (it == lp.iL[end]) + frc1[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) + 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 id1, id2 = lp.plidx[ipl2] Ush[b,1] = U[b,id1,r] Ush[b,2] = U[b,id2,r] sync_threads() + SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N) + bu1, ru1 = up((b, r), id1, lp) bu2, ru2 = up((b, r), id2, lp) bd, rd = up((bu1, ru1), id2, lp) @@ -368,16 +398,25 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::Spa if ru2 == r gt2 = Ush[bu2,1] else - gt2 = U[bu2,id1,ru2] + if SFBC && (it == lp.iL[end]) + gt2 = Ubnd + else + gt2 = U[bu2,id1,ru2] + end end l1 = gt1/gt2 l2 = Ush[b,2]\Ush[b,1] - 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) + if SFBC && (it == lp.iL[end]) + frc2[b,1,r] = projalg(Ush[b,1]*l1/Ush[b,2]) + 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 end diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 30657f9..9602e90 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -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 + ztw = ztwist(gp, lp) if abs(gp.c0-1) < 1.0E-10 @timeit "Wilson gauge action" 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 else @timeit "Improved gauge action" 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 @@ -38,9 +39,10 @@ end function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) + ztw = ztwist(gp, lp) @timeit "Plaquette measurement" 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 diff --git a/src/main/times.jl b/src/main/times.jl index b5abd0a..ce3069d 100644 --- a/src/main/times.jl +++ b/src/main/times.jl @@ -7,7 +7,8 @@ Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU") using LatticeGPU # 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) # Seed RNG @@ -21,6 +22,10 @@ ALG = SU3alg PREC = Float64 println("Precision: ", PREC) +# MD integrator +int = omf4(PREC, 0.1, 10) +println(int) + println("Allocating YM workspace") ymws = YMworkspace(GRP, PREC, lp) @@ -44,19 +49,17 @@ println("Initial Action: ") @time S = gauge_action(U, lp, gp, ymws) -dt = 0.1 -ns = 10 -HMC!(U,dt,1,lp, gp, ymws, noacc=true) +HMC!(U,int.eps,1,lp, gp, ymws, noacc=true) pl = Vector{Float64}() 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) push!(pl, plaquette(U,lp, gp, ymws)) println("# Plaquette: ", pl[end], "\n") - @time x, y = sfcoupling(U,lp,gp,ymws) - println("SF coupling: ", x, " ", y) +# @time x, y = sfcoupling(U,lp,gp,ymws) +# println("SF coupling: ", x, " ", y) end 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: ") @time wfl_rk3(U, 100, 0.01, gp, lp, ymws) eoft = Eoft_plaq(U, gp, lp, ymws) -eoft = Eoft_clover(U, lp, ymws) -qtop = Qtop(U, lp, ymws) +eoft = Eoft_clover(U, gp, lp, ymws) +qtop = Qtop(U, gp, lp, ymws) @time eoft = Eoft_plaq(U, gp, lp, ymws) println("Plaq: ", eoft) -@time eoft = Eoft_clover(U, lp, ymws) +@time eoft = Eoft_clover(U, gp, lp, ymws) println("Clov: ", eoft) -@time qtop = Qtop(U, lp, ymws) +@time qtop = Qtop(U, gp, lp, ymws) println("Qtop: ", qtop) println("Action: ", gauge_action(U, lp, gp, ymws)) @@ -88,13 +91,10 @@ println("Initial Action: ") @time S = gauge_action(U, lp, gp, ymws) -dt = 0.1 -ns = 10 -HMC!(U,dt,1,lp, gp, ymws, noacc=true) - +HMC!(U,int.eps,1,lp, gp, ymws, noacc=true) pl = Vector{Float64}() 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) push!(pl, plaquette(U,lp, gp, ymws)) println("# Plaquette: ", pl[end], "\n")