From ee1776c66518e3fd35d74a9e87a9b354f20d5bbd Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Sun, 24 Oct 2021 16:47:41 +0200 Subject: [PATCH] SF boundary conditions for Wilson action conserve energy --- src/YM/YM.jl | 35 +++++++++++++++++--- src/YM/YMact.jl | 81 +++++++++++++++++++++++++++++++++------------- src/YM/YMfields.jl | 27 +++++++++++++++- src/YM/YMflow.jl | 20 ++++++------ src/YM/YMhmc.jl | 8 ++--- src/main/times.jl | 14 ++++---- 6 files changed, 136 insertions(+), 49 deletions(-) diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 46d394c..9699d74 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -20,18 +20,43 @@ using ..MD import Base.show -struct GaugeParm{T} +struct GaugeParm{T,G} beta::T c0::T cG::NTuple{2,T} ng::Int64 + + Ubnd::G + + function GaugeParm{T}(::Type{G}, bt, c0, cG) where {T,G} + + function degree(::SU2{T}) where T + return 2 + end + function degree(::Type{SU3{T}}) where T <: AbstractFloat + return 3 + end + ng = degree(G) + + return new{T,G}(bt, c0, cG, ng, one(G)) + end + function GaugeParm{T}(::Type{G}, bt, c0) where {T,G} + + degree(::Type{SU2{T}}) where T <: AbstractFloat = 2 + degree(::Type{SU3{T}}) where T <: AbstractFloat = 3 + ng = degree(G) + + return new{T,G}(bt, c0, (0.0,0.0), ng, one(G)) + end end export GaugeParm -function Base.show(io::IO, gp::GaugeParm) +function Base.show(io::IO, gp::GaugeParm{T, G}) where {T,G} - println(io, "beta: ", gp.beta) - println(io, "c0: ", gp.c0) - println(io, "Ngauge: ", gp.ng) + 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) return nothing end diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 2907f27..7d50078 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -95,18 +95,22 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIO return nothing end -function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D} +function krnl_plaq!(plx, U::AbstractArray{T}, 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) + + 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)) - for id1 in 1:N-1 + for id1 in N:-1:1 bu1, ru1 = up((b, r), id1, lp) 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() @@ -114,15 +118,23 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,BC_PERIODIC,D}) if ru1 == r gt1 = Ush[bu1,2] else - gt1 = U[bu1,id2,ru1] + if SFBC && (it == lp.iL[end]) && (id1 == N) + gt1 = Ubnd + else + gt1 = U[bu1,id2,ru1] + end end if ru2 == r gt2 = Ush[bu2,1] else gt2 = U[bu2,id1,ru2] end - - S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + + if ITBND && SFBC && (id1 == N) + S += cG*tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + else + S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) + end end end @@ -132,9 +144,13 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,BC_PERIODIC,D}) return nothing end -function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D} +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} 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)) @@ -155,20 +171,39 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceP if ru1 == r gt1 = Ush[bu1,2] else - gt1 = U[bu1,id2,ru1] + if SFBC && (it == lp.iL[end]) && (id1 == N) + gt1 = Ubnd + else + gt1 = U[bu1,id2,ru1] + end end g1 = gt1/gt2 g2 = Ush[b,2]\Ush[b,1] - - X = projalg(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) + + if SFBC && (it == 1) && (id1 == N) + 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) + X = cG*projalg(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) + else + X = projalg(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) + end + end - + return nothing end @@ -308,31 +343,33 @@ end Computes the force deriving from the Wilson plaquette action, without the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` """ -function force_gauge(ymws::YMworkspace, U, c0, lp::SpaceParm) +function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm) if abs(c0-1) < 1.0E-10 @timeit "Wilson gauge force" begin - force_wilson_pln!(ymws.frc1, ymws.frc2, U, lp::SpaceParm) + force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, lp::SpaceParm) end else @timeit "Improved gauge force" begin - force_wilson_pln!(ymws.frc1, ymws.frc2, U, lp::SpaceParm, c0) + force_pln!(ymws.frc1, ymws.frc2, U, nothing, nothing, lp::SpaceParm, c0) end end return nothing end -force_wilson(ymws::YMworkspace, U, lp::SpaceParm) = force_gauge(ymws, U, 1, lp) +force_gauge(ymws::YMworkspace, U, c0, gp, lp) = force_gauge(ymws, U, c0, gp.cG[1], gp, lp) +force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp, lp) +force_wilson(ymws::YMworkspace, U, cG, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp.cG[1], gp, lp) -function force_wilson_pln!(frc1, ftmp, U, lp::SpaceParm, c0=1) +function force_pln!(frc1, ftmp, U, Ubnd, cG, 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,i,lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_wilson_pln!(frc1,ftmp,U, Ubnd, cG,i,lp) end end else diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl index df8094a..f5e043b 100644 --- a/src/YM/YMfields.jl +++ b/src/YM/YMfields.jl @@ -34,7 +34,7 @@ function randomize!(f, lp::SpaceParm, ymws::YMworkspace) return nothing end -function krnl_assign_SU3!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D} +function krnl_assign_SU3!(frc::AbstractArray{T}, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x for id in 1:lp.ndim @@ -45,6 +45,31 @@ function krnl_assign_SU3!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D return nothing end +function krnl_assign_SU3!(frc::AbstractArray{T}, m, 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) + + if ((B==BC_SF_AFWB)||(B==BC_SF_ORBI)) + if it == 1 + for id in 1:lp.ndim-1 + frc[b,id,r] = zero(T) + end + frc[b,N,r] = SU3alg(m[b,N,1,r], m[b,N,2,r], m[b,N,3,r], + m[b,N,4,r], m[b,N,5,r], m[b,N,6,r], + m[b,N,7,r], m[b,N,8,r]) + else + for id in 1:lp.ndim + frc[b,id,r] = SU3alg(m[b,id,1,r], m[b,id,2,r], m[b,id,3,r], + m[b,id,4,r], m[b,id,5,r], m[b,id,6,r], + m[b,id,7,r], m[b,id,8,r]) + end + end + end + + return nothing +end + function krnl_assign_SU2!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 5d4f959..85cac2c 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -61,11 +61,11 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S end -function flw_euler(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false) +function flw_euler(U, ns, eps, c0, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace; add_zth=false) @timeit "Integrating flow equations (Euler)" begin for i in 1:ns - force_gauge(ymws, U, c0, lp) + force_gauge(ymws, U, c0, 1, gp, lp) if add_zth add_zth_term(ymws::YMworkspace, U, lp) end @@ -76,12 +76,12 @@ function flw_euler(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=fal return nothing end -function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false) +function flw_rk3(U, ns, eps, c0, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace; add_zth=false) @timeit "Integrating flow equations (RK3)" begin for i in 1:ns e0 = eps/2 - force_gauge(ymws, U, c0, lp) + force_gauge(ymws, U, c0, 1, gp, lp) if add_zth add_zth_term(ymws::YMworkspace, U, lp) end @@ -90,7 +90,7 @@ function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false e0 = -34*eps/36 e1 = 16*eps/9 - force_gauge(ymws, U, c0, lp) + force_gauge(ymws, U, c0, 1, gp, lp) if add_zth add_zth_term(ymws::YMworkspace, U, lp) end @@ -98,7 +98,7 @@ function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false U .= expm.(U, ymws.mom) e1 = 6*eps/4 - force_gauge(ymws, U, c0, lp) + force_gauge(ymws, U, c0, 1, gp, lp) if add_zth add_zth_term(ymws::YMworkspace, U, lp) end @@ -113,10 +113,10 @@ end -wfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 1, lp, ymws) -zfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 5.0/3.0, lp, ymws, add_zth=true) -wfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 1, lp, ymws) -zfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 5.0/3.0, lp, ymws, add_zth=true) +wfl_euler(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 1, gp, lp, ymws) +zfl_euler(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 5.0/3.0, gp, lp, ymws, add_zth=true) +wfl_rk3(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 1, gp, lp, ymws) +zfl_rk3(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 5.0/3.0, gp, lp, ymws, add_zth=true) ## diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 751373d..e9216f4 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -20,7 +20,7 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) 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, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, gp.cG[1], lp) end end else @@ -40,7 +40,7 @@ function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) @timeit "Plaquette measurement" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, one(gp.cG[1]), lp) end end @@ -93,7 +93,7 @@ function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ym @timeit "MD evolution" begin ee = int.eps*gp.beta/gp.ng - force_gauge(ymws, U, gp.c0, lp) + force_gauge(ymws, U, gp.c0, gp, lp) mom .= mom .+ (int.r[1]*ee) .* ymws.frc1 for i in 1:int.ns k = 2 @@ -105,7 +105,7 @@ function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ym end k += off - force_gauge(ymws, U, gp.c0, lp) + force_gauge(ymws, U, gp.c0, gp, lp) if (i < int.ns) && (k == 1) mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1 else diff --git a/src/main/times.jl b/src/main/times.jl index 7076e95..a77a22a 100644 --- a/src/main/times.jl +++ b/src/main/times.jl @@ -7,7 +7,7 @@ Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU") using LatticeGPU # Set lattice/block size -lp = SpaceParm{4}((16,16,16,16), (4,4,4,4)) +lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_SF_AFWB) println("Space Parameters: ", lp) # Seed RNG @@ -36,7 +36,7 @@ println("Time to take the configuration to memory: ") # Set gauge parameters # FIRST SET: Wilson action/flow println("\n## WILSON ACTION/FLOW TIMES") -gp = GaugeParm{PREC}(6.0, 1.0, (0.0,0.0), 3) +gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0, (0.5,0.5)) println("Gauge Parameters: ", gp) @@ -56,11 +56,11 @@ for i in 1:4 println("# Plaquette: ", pl[end], "\n") end -wfl_rk3(U, 1, 0.01, lp, ymws) +wfl_rk3(U, 1, 0.01, gp, lp, ymws) println("Action: ", gauge_action(U, lp, gp, ymws)) println("Time for 100 steps of RK3 flow integrator: ") -@time wfl_rk3(U, 100, 0.01, lp, ymws) +@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) @@ -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}(6.0, 5.0/3.0, (0.0,0.0), 3) +gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0) println("Gauge Parameters: ", gp) println("Initial Action: ") @@ -97,11 +97,11 @@ for i in 1:4 println("# Plaquette: ", pl[end], "\n") end -zfl_rk3(U, 1, 0.01, lp, ymws) +zfl_rk3(U, 1, 0.01, gp, lp, ymws) println("Action: ", gauge_action(U, lp, gp, ymws)) println("Time for 100 steps of RK3 flow integrator: ") -@time zfl_rk3(U, 100, 0.01, lp, ymws) +@time zfl_rk3(U, 100, 0.01, gp, lp, ymws) println("Action: ", gauge_action(U, lp, gp, ymws)) println("## END improved action/flow measurements")