Working HMC for SF (w/wo improved actions).

Only Aoki-Frezzotti-Weisz Choice B available for improved actions

src/YM/YMact.jl

@@ -9,18 +9,20 @@
 ### created: Mon Jul 12 18:31:19 2021
 ###                               
 
-function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
+function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::T, cG, 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))
-    for id1 in 1:N-1
+    for id1 in N:-1:1
         bu1, ru1 = up((b, r), id1, lp)
+        SFBC  = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1==N) 
         Ush[b,1] = U[b,id1,r]
 
-        for id2 = id1+1:N
+        for id2 = 1:id1-1
             bu2, ru2 = up((b, r), id2, lp)
             Ush[b,2] = U[b,id2,r]
             sync_threads()
@@ -67,7 +69,11 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIO
             if r2 == r
                 gc = Ush[b2,2]
             else
-                gc = U[b2,id2,r2]
+                if SFBC && (it == lp.iL[end])
+                    gc = Ubnd
+                else
+                    gc = U[b2,id2,r2]
+                end
             end
             h3 = (ga*gb)/gc
             # END staples
@@ -80,12 +86,23 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIO
             if ru1 == r
                 ga = Ush[bu1,2]
             else
-                ga = U[bu1,id2,ru1]
+                if SFBC && (it == lp.iL[end])
+                    ga = Ubnd
+                else
+                    ga = U[bu1,id2,ru1]
+                end
             end
 
             g2 = Ush[b,2]\Ush[b,1]
             
-            S += c0*tr(g2*ga/gb) + c1*( tr(g2*h2/gb) + tr(g2*ga/h3))
+            if (it == lp.iL[end]) && SFBC
+                S += cG*(c0*tr(g2*ga/gb) + (3*c1/2)*tr(g2*ga/h3))
+            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))
+            end
+
         end
     end
 
@@ -100,9 +117,6 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
     it = point_time((b, r), lp)
 
-    ITBND = (it == 1) || (it == lp.iL[end])
-    SFBC  = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) )
-
     Ush = @cuStaticSharedMem(T, (D,2))
     
     S = zero(eltype(plx))
@@ -110,6 +124,9 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
         bu1, ru1 = up((b, r), id1, lp)
         Ush[b,1] = U[b,id1,r]
 
+        SFBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) &&
+                  ( (it == 1) || (it == lp.iL[end])) ) && (id1 == N) 
+
         for id2 = 1:id1-1
             bu2, ru2 = up((b, r), id2, lp)
             Ush[b,2] = U[b,id2,r]
@@ -118,7 +135,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
             if ru1 == r
                 gt1 = Ush[bu1,2]
             else
-                if SFBC && (it == lp.iL[end]) && (id1 == N)
+                if SFBND && (it == lp.iL[end])
                     gt1 = Ubnd
                 else
                     gt1 = U[bu1,id2,ru1]
@@ -130,7 +147,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd::T, cG, lp::SpaceParm{N,M,B,D
                 gt2 = U[bu2,id1,ru2]
             end
 
-            if ITBND && SFBC && (id1 == N)
+            if SFBND
                 S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
             else
                 S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
@@ -149,9 +166,6 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
     it = point_time((b, r), lp)
 
-    ITBND = (it == 1) || (it == lp.iL[end])
-    SFBC  = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) )
-
     Ush = @cuStaticSharedMem(T, (D,2))
     
     @inbounds begin
@@ -159,6 +173,8 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
         bu1, ru1 = up((b, r), id1, lp)
         bu2, ru2 = up((b, r), id2, lp)
 
+        SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N)
+
         Ush[b,1] = U[b,id1,r]
         Ush[b,2] = U[b,id2,r]
         sync_threads()
@@ -171,7 +187,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
         if ru1 == r
             gt1 = Ush[bu1,2]
         else
-            if SFBC && (it == lp.iL[end]) && (id1 == N)
+            if SFBC && (it == lp.iL[end])
                 gt1 = Ubnd
             else
                 gt1 = U[bu1,id2,ru1]
@@ -181,13 +197,13 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
         g1 = gt1/gt2
         g2 = Ush[b,2]\Ush[b,1]
 
-        if SFBC && (it == 1) && (id1 == N)
+        if SFBC && (it == 1)
             X = cG*projalg(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)
-        elseif SFBC && (it == lp.iL[end]) && (id1 == N)
+        elseif SFBC && (it == lp.iL[end])
             X = cG*projalg(Ush[b,1]*g1/Ush[b,2])
             
             frc1[b  ,id1, r ] -= X
@@ -207,9 +223,10 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd::T, cG, ip
     return nothing
 end
 
-function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
+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}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
+    it = point_time((b, r), lp)
 
     Ush = @cuStaticSharedMem(T, (D,2))
     
@@ -218,6 +235,8 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::
         bu1, ru1 = up((b, r), id1, lp)
         bu2, ru2 = up((b, r), id2, lp)
 
+        SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && (id1 == N)
+
         Ush[b,1] = U[b,id1,r]
         Ush[b,2] = U[b,id2,r]
         sync_threads()
@@ -236,7 +255,11 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::
         if r2 == r
             gc = Ush[b2,2]
         else
-            gc = U[b2,id2,r2]
+            if SFBC && (it == lp.iL[end])
+                gc = Ubnd
+            else
+                gc = U[b2,id2,r2]
+            end
         end
         h1 = (ga\gb)*gc
         
@@ -282,7 +305,11 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::
         if r2 == r
             gc = Ush[b2,2]
         else
-            gc = U[b2,id2,r2]
+            if SFBC && (it == lp.iL[end])
+                gc = Ubnd
+            else
+                gc = U[b2,id2,r2]
+            end
         end
         h3 = (ga*gb)/gc
             
@@ -313,25 +340,52 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::
         if ru1 == r
             ga = Ush[bu1,2]
         else
-            ga = U[bu1,id2,ru1]
+            if SFBC && (it == lp.iL[end])
+                ga = Ubnd
+            else
+                ga = U[bu1,id2,ru1]
+            end
         end
         
         g1 = ga/gb
         g2 = Ush[b,2]\Ush[b,1]
 
-        X = c0*projalg(Ush[b,1]*g1/Ush[b,2]) + c1 * (projalg(Ush[b,1]*h2/(Ush[b,2]*gb)) +
+        if SFBC && (it == 1)
-                                                     projalg(Ush[b,1]*ga/(Ush[b,2]*h3)))
+            X = (cG*c0)*projalg(Ush[b,1]*g1/Ush[b,2]) + c1*projalg(Ush[b,1]*h2/(Ush[b,2]*gb)) +
+                (3*c1*cG/2)*projalg(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
-        frc1[b,id2,r] += X + c1*projalg(h1*g1/Ush[b,2]) 
             
-        frc2[bu1,id2,ru1] -= c0*projalg(g1*g2) + c1*( projalg((ga/h3)*g2) +
+            frc2[bu1,id2,ru1] -= (cG*c0)*projalg(g1*g2) + (3*c1*cG/2)*projalg((ga/h3)*g2) +
-                                                      projalg((g1/h4)*Ush[b,1]) +
+                (3*c1*cG/2)*projalg((g1/Ush[b,2])*h1)
-                                                      projalg((g1/Ush[b,2])*h1) )
+            
+            frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + (3*c1*cG/2) * projalg((Ush[b,2]\h1)*g1) +
+                c1*projalg(g2*h2/gb) 
+        elseif SFBC && (it == lp.iL[end])
+            X = (cG*c0)*projalg(Ush[b,1]*g1/Ush[b,2]) +
+                (3*c1*cG/2) * (projalg(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 + (3*c1*cG/2)*projalg(h1*g1/Ush[b,2]) 
+            
+            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)))
+            
+            frc1[b,id1,r] -= X + c1*projalg(Ush[b,1]*g1/h4) 
+            frc1[b,id2,r] += X + c1*projalg(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[bu2,id1,ru2] += c0*projalg(g2*g1) + c1*( projalg((Ush[b,2]\h1)*g1) +
+                                                          projalg(g2*h2/gb) +
+                                                          projalg(h4\Ush[b,1]*g1) )
+        end
             
-        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) )
     end
         
     return nothing
@@ -351,7 +405,7 @@ function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm)
         end
     else
         @timeit "Improved gauge force" begin
-            force_pln!(ymws.frc1, ymws.frc2, U, nothing, nothing, lp::SpaceParm, c0)
+            force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, lp::SpaceParm, c0)
         end
     end
     return nothing
@@ -375,7 +429,7 @@ function force_pln!(frc1, ftmp, U, Ubnd, cG, lp::SpaceParm, c0=1)
     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,i,lp)
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_impr_pln!(frc1,ftmp,U,c0,(1-c0)/8,Ubnd, cG, i,lp)
             end
         end
     end

src/YM/YMhmc.jl

@@ -26,7 +26,7 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T})
     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, 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], lp)
             end
         end
     end

src/main/times.jl

@@ -78,7 +78,7 @@ println("## END Wilson action/flow measurements")
 # Set gauge parameters
 # SECOND SET: Improved action/flow
 println("\n## IMPROVED ACTION/FLOW TIMES")
-gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0)
+gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 5/6, (0.5,0.5))
 println("Gauge  Parameters: ", gp)
 
 println("Initial Action: ")