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

Implemented corrected TWBC for Wilson action/flow

Browse source
This commit is contained in:
Alberto Ramos 2023-01-28 10:56:05 +01:00
parent 4944740ad2
commit 9ba55b10f8
2 changed files with 27 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

23
src/YM/YMact.jl
View file

@ -84,14 +84,17 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, cG, ztw, lp::SpaceParm{N,M,B
@inbounds begin @inbounds begin
b = Int64(CUDA.threadIdx().x) b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x) r = Int64(CUDA.blockIdx().x)
TWP = (b==1) && (r==1)
it = point_time((b, r), lp) it = point_time((b, r), lp)
IBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) &&
( (it == 1) || (it == lp.iL[end])) )
S = zero(eltype(plx)) S = zero(eltype(plx))
ipl = 0 ipl = 0
for id1 in N:-1:1 for id1 in N:-1:1
bu1, ru1 = up((b, r), id1, lp) bu1, ru1 = up((b, r), id1, lp)
SFBND = ( ( (B == BC_SF_AFWB) || (B == BC_SF_ORBI) ) && SFBND = IBND && (id1 == N)
( (it == 1) || (it == lp.iL[end])) ) && (id1 == N)
for id2 = 1:id1-1 for id2 = 1:id1-1
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
@ -106,7 +109,11 @@ function krnl_plaq!(plx, U::AbstractArray{T}, Ubnd, cG, ztw, lp::SpaceParm{N,M,B
if SFBND if SFBND
S += cG*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) S += cG*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2]))
else else
if TWP
S += ztw[ipl]*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2])) S += ztw[ipl]*tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2]))
else
S += tr(U[b,id1,r]*gt1 / (U[b,id2,r]*U[bu2,id1,ru2]))
end
end end
end end
end end
@ -122,6 +129,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw,
b = Int64(CUDA.threadIdx().x) b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x) r = Int64(CUDA.blockIdx().x)
TWP = (b==1) && (r==1)
it = point_time((b,r), lp) it = point_time((b,r), lp)
@inbounds begin @inbounds begin
@ -153,12 +161,17 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw,
frc1[b ,id2, r ] += X frc1[b ,id2, r ] += X
frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1) frc2[bu2,id1,ru2] += cG*projalg(ztw,g2*g1)
else else
if TWP
X = projalg(ztw,U[b,id1,r]*g1/U[b,id2,r]) X = projalg(ztw,U[b,id1,r]*g1/U[b,id2,r])
frc1[b ,id1, r ] -= X
frc1[b ,id2, r ] += X
frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2) frc2[bu1,id2,ru1] -= projalg(ztw,g1*g2)
frc2[bu2,id1,ru2] += projalg(ztw,g2*g1) frc2[bu2,id1,ru2] += projalg(ztw,g2*g1)
else
X = projalg(U[b,id1,r]*g1/U[b,id2,r])
frc2[bu1,id2,ru1] -= projalg(g1*g2)
frc2[bu2,id1,ru2] += projalg(g2*g1)
end
frc1[b ,id1, r ] -= X
frc1[b ,id2, r ] += X
end end
end end

10
src/main/times.jl
View file

@ -7,7 +7,7 @@ Pkg.activate("../..")
using LatticeGPU using LatticeGPU
# Set lattice/block size # Set lattice/block size
ntwist = (0,0,0,0,0,0) ntwist = (1,0,0,0,0,0)
lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist) lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist)
println("Space Parameters: ", lp) println("Space Parameters: ", lp)
@ -41,10 +41,10 @@ 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}(GRP{PREC}, 6.0, 1.66666666) gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0)
println("Gauge Parameters: ", gp) println("Gauge Parameters: ", gp)
flwint = zfl_rk3(PREC, 0.005, 1.0E-6) flwint = wfl_rk3(PREC, 0.005, 1.0E-6)
println(flwint) println(flwint)
println("Initial Action: ") println("Initial Action: ")
@ -54,7 +54,7 @@ println("Initial Action: ")
HMC!(U,int.eps,1,lp, gp, ymws, noacc=true) HMC!(U,int.eps,1,lp, gp, ymws, noacc=true)
pl = Vector{Float64}() pl = Vector{Float64}()
for i in 1:100 for i in 1:4
@time dh, acc = HMC!(U,int,lp, gp, ymws) @time dh, acc = HMC!(U,int,lp, gp, ymws)
println("# HMC: ", acc, " ", dh) println("# HMC: ", acc, " ", dh)
push!(pl, plaquette(U,lp, gp, ymws)) push!(pl, plaquette(U,lp, gp, ymws))
@ -69,7 +69,7 @@ println("Action: ", gauge_action(U, lp, gp, ymws))
flw(U, flwint, 1, gp, lp, ymws) flw(U, flwint, 1, gp, lp, ymws)
println("Time for 100 steps of RK3 flow integrator: ") println("Time for 100 steps of RK3 flow integrator: ")
@time flw(U, flwint, 400, gp, lp, ymws) @time flw(U, flwint, 100, gp, lp, ymws)
eoft = Eoft_plaq(U, gp, lp, ymws) eoft = Eoft_plaq(U, gp, lp, ymws)
eoft = Eoft_clover(U, gp, lp, ymws) eoft = Eoft_clover(U, gp, lp, ymws)
qtop = Qtop(U, gp, lp, ymws) qtop = Qtop(U, gp, lp, ymws)