Added support for SF boundary field

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

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
+        degree(::Type{SU2{T}}) where T <: AbstractFloat = 2
-            return 2
+        degree(::Type{SU3{T}}) where T <: AbstractFloat = 3
-        end
-        function degree(::Type{SU3{T}}) where T <: AbstractFloat
-            return 3
-        end
         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

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

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)
@@ -468,18 +468,18 @@ function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, ipl1, ipl2
     bd, rd   = up((bu1, ru1), id2, lp)
     if ru1 == r
         gt1 = Ush[bu1,2]
+    else
+        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
-    end
     
     l1 = gt1/gt2
     l2 = Ush[b,2]\Ush[b,1]

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

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))