lattice/latticegpu.jl
1
0
Fork 0
mirror of https://igit.ific.uv.es/alramos/latticegpu.jl.git synced 2025-05-14 19:23:42 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

SF boundary conditions for Wilson action conserve energy

Browse source
This commit is contained in:
Alberto Ramos 2021-10-24 16:47:41 +02:00
parent 74e22502e3
commit ee1776c665
6 changed files with 136 additions and 49 deletions
Download patch file Download diff file Expand all files Collapse all files

35
src/YM/YM.jl
View file

@ -20,18 +20,43 @@ using ..MD
import Base.show import Base.show
struct GaugeParm{T} struct GaugeParm{T,G}
beta::T beta::T
c0::T c0::T
cG::NTuple{2,T} cG::NTuple{2,T}
ng::Int64 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 end
export GaugeParm 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, "Group: ", G)
println(io, "c0: ", gp.c0) println(io, " - beta: ", gp.beta)
println(io, "Ngauge: ", gp.ng) println(io, " - c0: ", gp.c0)
println(io, " - cG: ", gp.cG)
println(io, " - Boundary link: ", gp.Ubnd)
return nothing return nothing
end end

57
src/YM/YMact.jl
View file

@ -95,36 +95,48 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIO
return nothing return nothing
end 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 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)) Ush = @cuStaticSharedMem(T, (D,2))
S = zero(eltype(plx)) S = zero(eltype(plx))
for id1 in 1:N-1 for id1 in N:-1:1
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
Ush[b,1] = U[b,id1,r] Ush[b,1] = U[b,id1,r]
for id2 = id1+1:N for id2 = 1:id1-1
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
Ush[b,2] = U[b,id2,r] Ush[b,2] = U[b,id2,r]
sync_threads() sync_threads()
if ru1 == r if ru1 == r
gt1 = Ush[bu1,2] gt1 = Ush[bu1,2]
else
if SFBC && (it == lp.iL[end]) && (id1 == N)
gt1 = Ubnd
else else
gt1 = U[bu1,id2,ru1] gt1 = U[bu1,id2,ru1]
end end
end
if ru2 == r if ru2 == r
gt2 = Ush[bu2,1] gt2 = Ush[bu2,1]
else else
gt2 = U[bu2,id1,ru2] gt2 = U[bu2,id1,ru2]
end end
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)) S += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
end end
end end
end
I = point_coord((b,r), lp) I = point_coord((b,r), lp)
plx[I] = S plx[I] = S
@ -132,9 +144,13 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,BC_PERIODIC,D})
return nothing return nothing
end 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 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)) Ush = @cuStaticSharedMem(T, (D,2))
@ -154,13 +170,30 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceP
end end
if ru1 == r if ru1 == r
gt1 = Ush[bu1,2] gt1 = Ush[bu1,2]
else
if SFBC && (it == lp.iL[end]) && (id1 == N)
gt1 = Ubnd
else else
gt1 = U[bu1,id2,ru1] gt1 = U[bu1,id2,ru1]
end end
end
g1 = gt1/gt2 g1 = gt1/gt2
g2 = Ush[b,2]\Ush[b,1] g2 = Ush[b,2]\Ush[b,1]
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]) X = projalg(Ush[b,1]*g1/Ush[b,2])
frc1[b ,id1, r ] -= X frc1[b ,id1, r ] -= X
@ -169,6 +202,8 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceP
frc2[bu2,id1,ru2] += projalg(g2*g1) frc2[bu2,id1,ru2] += projalg(g2*g1)
end end
end
return nothing return nothing
end end
@ -308,31 +343,33 @@ end
Computes the force deriving from the Wilson plaquette action, without Computes the force deriving from the Wilson plaquette action, without
the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` 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 if abs(c0-1) < 1.0E-10
@timeit "Wilson gauge force" begin @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 end
else else
@timeit "Improved gauge force" begin @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
end end
return nothing return nothing
end 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!(frc1, zero(eltype(frc1)))
fill!(ftmp, zero(eltype(ftmp))) fill!(ftmp, zero(eltype(ftmp)))
if c0 == 1 if c0 == 1
for i in 1:lp.npls for i in 1:lp.npls
CUDA.@sync begin 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
end end
else else

27
src/YM/YMfields.jl
View file

@ -34,7 +34,7 @@ function randomize!(f, lp::SpaceParm, ymws::YMworkspace)
return nothing return nothing
end 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 b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
for id in 1:lp.ndim 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 return nothing
end 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} 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 b, r = CUDA.threadIdx().x, CUDA.blockIdx().x

20
src/YM/YMflow.jl
View file

@ -61,11 +61,11 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S
end 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 @timeit "Integrating flow equations (Euler)" begin
for i in 1:ns for i in 1:ns
force_gauge(ymws, U, c0, lp) force_gauge(ymws, U, c0, 1, gp, lp)
if add_zth if add_zth
add_zth_term(ymws::YMworkspace, U, lp) add_zth_term(ymws::YMworkspace, U, lp)
end end
@ -76,12 +76,12 @@ function flw_euler(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=fal
return nothing return nothing
end 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 @timeit "Integrating flow equations (RK3)" begin
for i in 1:ns for i in 1:ns
e0 = eps/2 e0 = eps/2
force_gauge(ymws, U, c0, lp) force_gauge(ymws, U, c0, 1, gp, lp)
if add_zth if add_zth
add_zth_term(ymws::YMworkspace, U, lp) add_zth_term(ymws::YMworkspace, U, lp)
end end
@ -90,7 +90,7 @@ function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false
e0 = -34*eps/36 e0 = -34*eps/36
e1 = 16*eps/9 e1 = 16*eps/9
force_gauge(ymws, U, c0, lp) force_gauge(ymws, U, c0, 1, gp, lp)
if add_zth if add_zth
add_zth_term(ymws::YMworkspace, U, lp) add_zth_term(ymws::YMworkspace, U, lp)
end end
@ -98,7 +98,7 @@ function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false
U .= expm.(U, ymws.mom) U .= expm.(U, ymws.mom)
e1 = 6*eps/4 e1 = 6*eps/4
force_gauge(ymws, U, c0, lp) force_gauge(ymws, U, c0, 1, gp, lp)
if add_zth if add_zth
add_zth_term(ymws::YMworkspace, U, lp) add_zth_term(ymws::YMworkspace, U, lp)
end end
@ -113,10 +113,10 @@ end
wfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 1, lp, ymws) 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, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 5.0/3.0, lp, ymws, add_zth=true) 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, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 1, lp, ymws) 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, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 5.0/3.0, lp, ymws, add_zth=true) 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)
## ##

8
src/YM/YMhmc.jl
View file

@ -20,7 +20,7 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T})
if abs(gp.c0-1) < 1.0E-10 if abs(gp.c0-1) < 1.0E-10
@timeit "Wilson gauge action" begin @timeit "Wilson gauge action" begin
CUDA.@sync 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
end end
else else
@ -40,7 +40,7 @@ function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace)
@timeit "Plaquette measurement" begin @timeit "Plaquette measurement" begin
CUDA.@sync 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
end end
@ -93,7 +93,7 @@ function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ym
@timeit "MD evolution" begin @timeit "MD evolution" begin
ee = int.eps*gp.beta/gp.ng 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 mom .= mom .+ (int.r[1]*ee) .* ymws.frc1
for i in 1:int.ns for i in 1:int.ns
k = 2 k = 2
@ -105,7 +105,7 @@ function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ym
end end
k += off k += off
force_gauge(ymws, U, gp.c0, lp) force_gauge(ymws, U, gp.c0, gp, lp)
if (i < int.ns) && (k == 1) if (i < int.ns) && (k == 1)
mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1 mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1
else else

14
src/main/times.jl
View file

@ -7,7 +7,7 @@ Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
using LatticeGPU using LatticeGPU
# Set lattice/block size # 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) println("Space Parameters: ", lp)
# Seed RNG # Seed RNG
@ -36,7 +36,7 @@ println("Time to take the configuration to memory: ")
# Set gauge parameters # Set gauge parameters
# FIRST SET: Wilson action/flow # FIRST SET: Wilson action/flow
println("\n## WILSON ACTION/FLOW TIMES") 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) println("Gauge Parameters: ", gp)
@ -56,11 +56,11 @@ for i in 1:4
println("# Plaquette: ", pl[end], "\n") println("# Plaquette: ", pl[end], "\n")
end 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("Action: ", gauge_action(U, lp, gp, ymws))
println("Time for 100 steps of RK3 flow integrator: ") 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_plaq(U, gp, lp, ymws)
eoft = Eoft_clover(U, lp, ymws) eoft = Eoft_clover(U, lp, ymws)
qtop = Qtop(U, lp, ymws) qtop = Qtop(U, lp, ymws)
@ -78,7 +78,7 @@ println("## END Wilson action/flow measurements")
# Set gauge parameters # Set gauge parameters
# SECOND SET: Improved action/flow # SECOND SET: Improved action/flow
println("\n## IMPROVED ACTION/FLOW TIMES") 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("Gauge Parameters: ", gp)
println("Initial Action: ") println("Initial Action: ")
@ -97,11 +97,11 @@ for i in 1:4
println("# Plaquette: ", pl[end], "\n") println("# Plaquette: ", pl[end], "\n")
end 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("Action: ", gauge_action(U, lp, gp, ymws))
println("Time for 100 steps of RK3 flow integrator: ") 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("Action: ", gauge_action(U, lp, gp, ymws))
println("## END improved action/flow measurements") println("## END improved action/flow measurements")