diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index ad2229d..17d6277 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -44,6 +44,6 @@ export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export Eoft_clover, Eoft_plaq, Qtop export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 export flw, flw_adapt -export sfcoupling +export sfcoupling, bndfield, setbndfield end # module diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 8288c0b..3b2d10a 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -20,25 +20,22 @@ using ..MD import Base.show -struct GaugeParm{T,G} +struct GaugeParm{T,G,N} beta::T c0::T cG::NTuple{2,T} ng::Int64 - Ubnd::G + Ubnd::NTuple{N, G} - function GaugeParm{T}(::Type{G}, bt, c0, cG) where {T,G} + function GaugeParm{T}(::Type{G}, bt, c0, cG, phi, iL) where {T,G} - function degree(::Type{SU2{T}}) where T - return 2 - end - function degree(::Type{SU3{T}}) where T <: AbstractFloat - return 3 - end + degree(::Type{SU2{T}}) where T <: AbstractFloat = 2 + degree(::Type{SU3{T}}) where T <: AbstractFloat = 3 ng = degree(G) + nsd = length(iL) - return new{T,G}(bt, c0, cG, ng, one(G)) + return new{T,G,nsd}(bt, c0, cG, ng, ntuple(id->bndfield(phi[1], phi[2], iL[id]), nsd)) end function GaugeParm{T}(::Type{G}, bt, c0) where {T,G} @@ -46,17 +43,21 @@ struct GaugeParm{T,G} degree(::Type{SU3{T}}) where T <: AbstractFloat = 3 ng = degree(G) - return new{T,G}(bt, c0, (0.0,0.0), ng, one(G)) + return new{T,G,0}(bt, c0, (0.0,0.0), ng, ()) end end export GaugeParm -function Base.show(io::IO, gp::GaugeParm{T, G}) where {T,G} +function Base.show(io::IO, gp::GaugeParm{T, G, N}) where {T,G,N} println(io, "Group: ", G) println(io, " - beta: ", gp.beta) println(io, " - c0: ", gp.c0) println(io, " - cG: ", gp.cG) - println(io, " - Boundary link: ", gp.Ubnd) + if (N > 0) + for i in 1:N + println(io, " - Boundary link: ", gp.Ubnd[i]) + end + end return nothing end @@ -142,6 +143,6 @@ export Eoft_clover, Eoft_plaq, Qtop export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 include("YMsf.jl") -export sfcoupling +export sfcoupling, bndfield, setbndfield end diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index be22107..9bdb835 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -9,7 +9,7 @@ ### created: Mon Jul 12 18:31:19 2021 ### -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} +function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::NTuple{NB,T}, cG, ztw, lp::SpaceParm{N,M,B,D}) where {T,NB,N,M,B,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x it = point_time((b, r), lp) @@ -72,7 +72,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ztw, lp::Spac gc = Ush[b2,2] else if SFBC && (it == lp.iL[end]) - gc = Ubnd + gc = Ubnd[id2] else gc = U[b2,id2,r2] end @@ -89,7 +89,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ztw, lp::Spac ga = Ush[bu1,2] else if SFBC && (it == lp.iL[end]) - ga = Ubnd + ga = Ubnd[id2] else ga = U[bu1,id2,ru1] end @@ -115,7 +115,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, ztw, lp::Spac return nothing end -function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, ztw, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} +function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, 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) @@ -141,7 +141,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, ztw, lp::SpaceParm{N, gt1 = Ush[bu1,2] else if SFBND && (it == lp.iL[end]) - gt1 = Ubnd + gt1 = Ubnd[id2] else gt1 = U[bu1,id2,ru1] end @@ -166,7 +166,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, ztw, lp::SpaceParm{N, return nothing end -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} +function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, 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) @@ -193,7 +193,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, zt gt1 = Ush[bu1,2] else if SFBC && (it == lp.iL[end]) - gt1 = Ubnd + gt1 = Ubnd[id2] else gt1 = U[bu1,id2,ru1] end @@ -227,7 +227,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, zt return nothing end -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} +function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, 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) @@ -260,7 +260,7 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T, gc = Ush[b2,2] else if SFBC && (it == lp.iL[end]) - gc = Ubnd + gc = Ubnd[id2] else gc = U[b2,id2,r2] end @@ -310,7 +310,7 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T, gc = Ush[b2,2] else if SFBC && (it == lp.iL[end]) - gc = Ubnd + gc = Ubnd[id2] else gc = U[b2,id2,r2] end @@ -345,7 +345,7 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd::T, ga = Ush[bu1,2] else if SFBC && (it == lp.iL[end]) - ga = Ubnd + ga = Ubnd[id2] else ga = U[bu1,id2,ru1] end diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 3235677..62eec71 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -259,7 +259,7 @@ 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}, Ubnd::T, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} +function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x @@ -270,7 +270,7 @@ function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd::T, ztw, ipl, lp::SpacePa bu2, ru2 = up((b, r), id2, lp) if SFBC && (ru1 != r) - gt = Ubnd + gt = Ubnd[id2] else gt = U[bu1,id2,ru1] end @@ -407,7 +407,7 @@ function krnl_add_qd!(rm, op, frc1, frc2, U, lp::SpaceParm{4,M,B,D}) where {M,B, return nothing end -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} +function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, Ubnd, 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) @@ -469,16 +469,16 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, ipl1, ipl2 if ru1 == r gt1 = Ush[bu1,2] else - gt1 = U[bu1,id2,ru1] + if SFBC && (it == lp.iL[end]) + gt2 = Ubnd[id2] + else + gt1 = U[bu1,id2,ru1] + end end if ru2 == r gt2 = Ush[bu2,1] else - if SFBC && (it == lp.iL[end]) - gt2 = Ubnd - else - gt2 = U[bu2,id1,ru2] - end + gt2 = U[bu2,id1,ru2] end l1 = gt1/gt2 diff --git a/src/YM/YMsf.jl b/src/YM/YMsf.jl index faa6f68..2bebb62 100644 --- a/src/YM/YMsf.jl +++ b/src/YM/YMsf.jl @@ -43,7 +43,7 @@ function sfcoupling(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) return dsdeta, ddnu end -function krnl_sfcoupling!(rm, U::AbstractArray{T}, Ubnd::T, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} +function krnl_sfcoupling!(rm, U::AbstractArray{T}, Ubnd, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x I = point_coord((b,r), lp) @@ -68,7 +68,7 @@ function krnl_sfcoupling!(rm, U::AbstractArray{T}, Ubnd::T, lp::SpaceParm{N,M,B, for id in 1:N-1 bu, ru = up((b,r), id, lp) - X = projalg(Ubnd/(U[b,id,r]*U[bu,N,ru])*U[b,N,r]) + X = projalg(Ubnd[id]/(U[b,id,r]*U[bu,N,ru])*U[b,N,r]) rm[I] -= 3*X.t7 + SR3 * X.t8 rm[ID] += 2*X.t7 - SR3x2 * X.t8 end @@ -76,3 +76,40 @@ function krnl_sfcoupling!(rm, U::AbstractArray{T}, Ubnd::T, lp::SpaceParm{N,M,B, return nothing end + + +@inline function bndfield(phi1::T, phi2::T, iL) where T <: AbstractFloat + + SR3::T = 1.73205080756887729352744634151 + + zt = zero(T) + X = SU3alg{T}(zt,zt,zt,zt,zt,zt,(phi1-phi2)/iL,SR3*(phi1+phi2)/iL) + + return exp(X) +end + + +function setbndfield(U, phi, lp::SpaceParm{N,M,B,D}) where {N,M,B,D} + + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_setbnd_it0!(U, phi[1,1], phi[1,2], lp) + end + + return nothing +end + +function krnl_setbnd_it0!(U, phi1, phi2, lp::SpaceParm) + + b, r = CUDA.threadIdx().x, CUDA.blockIdx().x + it = point_time((b,r), lp) + + SFBC = (B == BC_SF_AFWB) || (B == BC_SF_ORBI) + + if (it == 0) && SFBC + for id in 1:N-1 + U[b,id,r] = bndfield(phi1,phi2,lp.iL[id]) + end + end + + return nothing +end diff --git a/tests/test_su3.jl b/tests/test_su3.jl index 26aeb4a..47972c5 100644 --- a/tests/test_su3.jl +++ b/tests/test_su3.jl @@ -134,3 +134,12 @@ ba = rand(SU3alg{T}) ga = exp(ba) println("Matrix: ", alg2mat(ba)) println("Exp: ", ga) + + +println("## Final tests: ") +g1 = exp(SU3alg{T}(6.23, 1.23, -0.34, 2.34, -0.23, 0.23, -8.34, 8.34)) +println(g1) +println(isgroup(g1)) +g1 = unitarize(g1) +println(g1) +println(isgroup(g1))