From 5336189ad5ee19711d960cacd7fb44bcdd41f83a Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Mon, 25 Oct 2021 15:34:07 +0200 Subject: [PATCH] Working HMC for SF (w/wo improved actions). Only Aoki-Frezzotti-Weisz Choice B available for improved actions --- src/YM/YMact.jl | 126 +++++++++++++++++++++++++++++++++------------- src/YM/YMhmc.jl | 2 +- src/main/times.jl | 2 +- 3 files changed, 92 insertions(+), 38 deletions(-) diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 7d50078..8672791 100644 --- a/src/YM/YMact.jl +++ b/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)) + - 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) ) + if SFBC && (it == 1) + 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 + + frc2[bu1,id2,ru1] -= (cG*c0)*projalg(g1*g2) + (3*c1*cG/2)*projalg((ga/h3)*g2) + + (3*c1*cG/2)*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 + 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 diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index e9216f4..39e263e 100644 --- a/src/YM/YMhmc.jl +++ b/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 diff --git a/src/main/times.jl b/src/main/times.jl index a77a22a..47e4cfe 100644 --- a/src/main/times.jl +++ b/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: ")