fixed the OBC branch in the Zeuthen term - now results match openQCD

This commit is contained in:
Nicolas Lang 2025-05-26 12:48:04 +02:00
parent 58685f63a2
commit 9ba07b506a
2 changed files with 27 additions and 5 deletions

View file

@ -165,7 +165,7 @@ include("YMfields.jl")
export randomize!, zero!, norm2
include("YMact.jl")
export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson
export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson, bnd_rescale_flw!
include("YMhmc.jl")
export gauge_action, hamiltonian, plaquette, HMC!, MD!

View file

@ -127,6 +127,22 @@ function add_zth_term(ymws::YMworkspace, U, lp)
return nothing
end
"""
function rescale_bnd(ymws::YMworkspace, lp)
Rescales the time-boundary of ymws.frc1 by a factor of 2.0 (needed for open BC)
"""
function rescale_bnd(ymws::YMworkspace, lp)
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz bnd_rescale_flw!(ymws.frc1,lp::SpaceParm)
end
return nothing
end
function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::SpaceParm{N,M,B,D}) where {TA,TG,N,M,B,D}
@inbounds begin
@ -155,12 +171,14 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S
end
end
if OBC
if (it > 1) && (it < lp.iL[end])
if (id < N)
frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
projalg(U[b,id,r]*X/U[b,id,r]))
elseif ((it == lp.iL[end]) || (it == 1)) && (id < N)
elseif (it > 1) && (it < (lp.iL[end]-1))
frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
projalg(U[b,id,r]*X/U[b,id,r]))
elseif (it == 1) || (it == (lp.iL[end]-1))
frc2[b,id,r] = frc[b,id,r]
end
else
frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
@ -171,6 +189,8 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S
return nothing
end
"""
function flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace)
@ -204,18 +224,20 @@ flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMwor
function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D}
@timeit "Integrating flow equations" begin
for i in 1:ns
force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
force_gauge(ymws, U, int.c0, 1, gp, lp)
if int.add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
rescale_bnd(ymws, lp)
ymws.mom .= ymws.frc1
U .= expm.(U, ymws.mom, 2*eps*int.r)
for k in 1:NI
force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
force_gauge(ymws, U, int.c0, 1, gp, lp)
if int.add_zth
add_zth_term(ymws::YMworkspace, U, lp)
end
rescale_bnd(ymws, lp)
ymws.mom .= int.e0[k].*ymws.mom .+ int.e1[k].*ymws.frc1
U .= expm.(U, ymws.mom, 2*eps)
end