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

Bugs corrected

Browse source
This commit is contained in:
Alberto Ramos 2021-09-25 15:22:47 +02:00
parent b347d33dbf
commit e8493693b0
6 changed files with 35 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/LatticeGPU.jl
View file

@ -31,5 +31,6 @@ include("YM/YM.jl")
using .YM using .YM
export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2 export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2
export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
export wfl_euler, wfl_rk3
end # module end # module

3
src/YM/YM.jl
View file

@ -102,4 +102,7 @@ export krnl_plaq!, force0_wilson!
include("YMhmc.jl") include("YMhmc.jl")
export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
include("YMflow.jl")
export wfl_euler, wfl_rk3
end end

26
src/YM/YMact.jl
View file

@ -9,7 +9,7 @@
### created: Mon Jul 12 18:31:19 2021 ### created: Mon Jul 12 18:31:19 2021
### ###
function krnl_plaq!(plx, U, lp::SpaceParm{N,M,D}) where {T,N,M,D} function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
@ -17,25 +17,29 @@ function krnl_plaq!(plx, U, lp::SpaceParm{N,M,D}) where {T,N,M,D}
plx[b,r] = zero(plx[b,r]) plx[b,r] = zero(plx[b,r])
for id1 in 1:N-1 for id1 in 1:N-1
if ru2 == r bu1, ru1 = up((b, r), id1, lp)
gt2 = Ush[bu2,1] Ush[b,1] = U[b,id1,r]
else
gt2 = U[bu2,id1,ru2]
end
for id2 = id1+1:N for id2 = id1+1:N
bu2, ru2 = up((b, r), id2, lp)
Ush[b,2] = U[b,id2,r]
sync_threads()
if ru1 == r if ru1 == r
gt1 = Ush[bu1,2] gt1 = Ush[bu1,2]
else else
gt1 = U[bu1,id2,ru1] gt1 = U[bu1,id2,ru1]
end end
sync_threads() if ru2 == r
gt2 = Ush[bu2,1]
bu1, ru1 = up((b, r), id1, lp) else
bu2, ru2 = up((b, r), id2, lp) gt2 = U[bu2,id1,ru2]
end
plx[b,r] += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) plx[b,r] += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2))
end
end end
return nothing return nothing
end end

8
src/YM/YMflow.jl
View file

@ -13,7 +13,7 @@ function wfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
for i in 1:ns for i in 1:ns
force_wilson(ymws, U, lp) force_wilson(ymws, U, lp)
U .= expm(U, ymws.frc1, 2*eps) U .= expm.(U, ymws.frc1, 2*eps)
end end
return nothing return nothing
@ -25,18 +25,18 @@ function wfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
c0 = eps/2 c0 = eps/2
force_wilson(ymws, U, lp) force_wilson(ymws, U, lp)
ymws.mom .= ymws.frc1 ymws.mom .= ymws.frc1
U .= expm(U, ymws.mom, c0) U .= expm.(U, ymws.mom, c0)
c0 = -34*eps/36 c0 = -34*eps/36
c1 = 16*eps/9 c1 = 16*eps/9
force_wilson(ymws, U, lp) force_wilson(ymws, U, lp)
ymws.mom .= c0.*ymws.mom .+ c1.*ymws.frc1 ymws.mom .= c0.*ymws.mom .+ c1.*ymws.frc1
U .= expm(U, ymws.mom) U .= expm.(U, ymws.mom)
c1 = 6*eps/4 c1 = 6*eps/4
force_wilson(ymws, U, lp) force_wilson(ymws, U, lp)
ymws.mom .= c1.*ymws.frc1 .- ymws.mom ymws.mom .= c1.*ymws.frc1 .- ymws.mom
U .= expm(U, ymws.mom) U .= expm.(U, ymws.mom)
end end
return nothing return nothing

2
src/YM/YMhmc.jl
View file

@ -15,7 +15,7 @@
Returns the value of the Wilson plaquette action for the configuration U. Returns the value of the Wilson plaquette action for the configuration U.
""" """
function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: A function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat
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, lp)

16
src/main/test.jl
View file

@ -12,7 +12,7 @@ println(CUDA.device())
GRP = SU3 GRP = SU3
ALG = SU3alg ALG = SU3alg
PREC = Float64 PREC = Float64
lp = SpaceParm{4}((64,64,64,64), (4,4,4,4)) lp = SpaceParm{4}((32,32,32,32), (4,4,4,4))
gp = GaugeParm{PREC}(6.0, (0.0,0.0), 3) gp = GaugeParm{PREC}(6.0, (0.0,0.0), 3)
println("Space Parameters: ", lp) println("Space Parameters: ", lp)
@ -28,6 +28,10 @@ println("Allocating gauge field")
U = field(GRP{PREC}, lp) U = field(GRP{PREC}, lp)
fill!(U, one(GRP{PREC})) fill!(U, one(GRP{PREC}))
println("Take to take the configuration to memory: ")
@time Ucpu = Array(U)
println("Allocating YM workspace") println("Allocating YM workspace")
ymws = YMworkspace(GRP, PREC, lp, save_mem = true) ymws = YMworkspace(GRP, PREC, lp, save_mem = true)
@ -60,10 +64,12 @@ for i in 1:ntot
println("# Plaquette: ", pl[end], "\n") println("# Plaquette: ", pl[end], "\n")
end end
io = open("foo.txt", "w") @time wfl_rk3(U, 1, 0.01, lp, ymws)
for i in 1:ntot
println(io, pl[i]) println("Action: ", gauge_action(U, lp, gp, ymws))
end println("Time for 100 steps of RK3 flow integrator: ")
@time wfl_rk3(U, 100, 0.01, lp, ymws)
println("Action: ", gauge_action(U, lp, gp, ymws))
println("END") println("END")