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

Added routines for Zeuthen flow

Browse source
This commit is contained in:
Alberto Ramos 2021-09-30 12:59:00 +02:00
parent de3ee1c7f5
commit 334ddf813a
2 changed files with 121 additions and 46 deletions
Download patch file Download diff file Expand all files Collapse all files

121
src/YM/YMflow.jl
View file

@ -9,69 +9,106 @@
### created: Sat Sep 25 08:37:14 2021
###
function wfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
"""
function add_zth_term(ymws::YMworkspace, U, lp)
Assuming that the gauge improved (LW) force is in ymws.frc1, this routine
adds the "Zeuthen term" and returns the full zeuthen force in ymws.frc1
"""
function add_zth_term(ymws::YMworkspace, U, lp)
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_zth!(ymws.frc1,ymws.frc2,U,lp)
end
ymws.frc1 .= (5/6).*ymws.frc1 .+ ymws.frc2
return nothing
end
function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::SpaceParm{N,M,D}) where {TA,TG,N,M,D}
b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
Ush = @cuStaticSharedMem(TG, D)
Ush = @cuStaticSharedMem(TA, D)
@inbounds for id in 1:N
Ush[b] = U[b,id,r]
Fsh[b] = Frc[b,id,r]
sync_threads()
bu, ru, bd, rd = updw((b,r), id, lp)
if ru == r
X = Fsh[bu]
else
X = frc[bu,id,ru]
end
if rd == r
Y = Fsh[bd]
Ud = Ush[bd]
else
Y = frc[bd,id,rd]
Ud = U[bd,id,rd]
end
frc2[b,id,r] = (1/6)*(projalg(Ud[b]\Y*Ud[b]) +
projalg(Ush[b]*X/Ush[b]))
end
return nothing
end
function flw_euler(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false)
for i in 1:ns
force_wilson(ymws, U, lp)
force_gauge(ymws, U, c0, lp)
if add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
U .= expm.(U, ymws.frc1, 2*eps)
end
return nothing
end
function wfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
for i in 1:ns
c0 = eps/2
force_wilson(ymws, U, lp)
ymws.mom .= ymws.frc1
U .= expm.(U, ymws.mom, c0)
c0 = -34*eps/36
c1 = 16*eps/9
force_wilson(ymws, U, lp)
ymws.mom .= c0.*ymws.mom .+ c1.*ymws.frc1
U .= expm.(U, ymws.mom)
c1 = 6*eps/4
force_wilson(ymws, U, lp)
ymws.mom .= c1.*ymws.frc1 .- ymws.mom
U .= expm.(U, ymws.mom)
end
return nothing
end
function zfl_euler(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
for i in 1:ns
force_gauge(ymws, U, 5.0/3.0, lp)
U .= expm.(U, ymws.frc1, 2*eps)
end
function flw_rk3(U, ns, eps, c0, lp::SpaceParm, ymws::YMworkspace; add_zth=false)
return nothing
end
function zfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace)
for i in 1:ns
c0 = eps/2
force_gauge(ymws, U, 5.0/3.0, lp)
force_gauge(ymws, U, c0, lp)
if add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
ymws.mom .= ymws.frc1
U .= expm.(U, ymws.mom, c0)
c0 = -34*eps/36
c1 = 16*eps/9
force_gauge(ymws, U, 5.0/3.0, lp)
force_gauge(ymws, U, c0, lp)
if add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
ymws.mom .= c0.*ymws.mom .+ c1.*ymws.frc1
U .= expm.(U, ymws.mom)
c1 = 6*eps/4
force_gauge(ymws, U, 5.0/3.0, lp)
force_gauge(ymws, U, c0, lp)
if add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
ymws.mom .= c1.*ymws.frc1 .- ymws.mom
U .= expm.(U, ymws.mom)
end
return nothing
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)