From 830163e90f8e70d7c83913f82d1340586e0c607a Mon Sep 17 00:00:00 2001 From: Nicolas Lang Date: Mon, 27 Nov 2023 13:19:26 +0100 Subject: [PATCH 1/5] added rng arguments to HMC function --- src/YM/YMfields.jl | 6 +++--- src/YM/YMhmc.jl | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl index c6e5433..d853c49 100644 --- a/src/YM/YMfields.jl +++ b/src/YM/YMfields.jl @@ -9,11 +9,11 @@ ### created: Thu Jul 15 15:16:47 2021 ### -function randomize!(f, lp::SpaceParm, ymws::YMworkspace) +function randomize!(f, lp::SpaceParm, ymws::YMworkspace; curng=CUDA.default_rng()) if ymws.ALG == SU2alg @timeit "Randomize SU(2) algebra field" begin - m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,3,lp.rsz) + m = Random.randn(curng, ymws.PRC, lp.bsz,lp.ndim,3,lp.rsz) CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU2!(f,m,lp) end @@ -23,7 +23,7 @@ function randomize!(f, lp::SpaceParm, ymws::YMworkspace) if ymws.ALG == SU3alg @timeit "Randomize SU(3) algebra field" begin - m = CUDA.randn(ymws.PRC, lp.bsz,lp.ndim,8,lp.rsz) + m = Random.randn(curng, ymws.PRC, lp.bsz,lp.ndim,8,lp.rsz) CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU3!(f,m,lp) end diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index c973875..81a93c0 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -58,13 +58,13 @@ function hamiltonian(mom, U, lp, gp, ymws) return K+V end -function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T +function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false, rng=Random.default_rng(), curng=CUDA.default_rng()) where T @timeit "HMC trayectory" begin ymws.U1 .= U - randomize!(ymws.mom, lp, ymws) + randomize!(ymws.mom, lp, ymws; curng) hini = hamiltonian(ymws.mom, U, lp, gp, ymws) MD!(ymws.mom, U, int, lp, gp, ymws) @@ -78,7 +78,7 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac end if (pacc < 1.0) - r = rand() + r = rand(rng) if (pacc < r) U .= ymws.U1 acc = false @@ -90,7 +90,7 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac end return dh, acc end -HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc) +HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false, rng=Random.default_rng(), curng=CUDA.default_rng()) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc, rng, curng) function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat} From 3256fe917c87a1247168555dd3b3be038d53f627 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Mon, 1 Jan 2024 21:23:14 +0100 Subject: [PATCH 2/5] Added md source to doc --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index 890b3f1..f7b4cf1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,7 @@ makedocs(sitename="LatticeGPU", modules=[LatticeGPU], doctest=true, "Yang-Mills" => "ym.md", "Gradient flow" => "flow.md", "Schrödinger Functional" => "sf.md", + "Spinors" => "spinors.md", "Dirac" => "dirac.md", "Solvers" => "solvers.md", "Input Output" => "io.md" From 094390c1eea6142799e508fa8c1dd59869726c2c Mon Sep 17 00:00:00 2001 From: Nicolas Lang Date: Fri, 2 May 2025 16:30:19 +0200 Subject: [PATCH 3/5] factor of 1/2 on spatial bnd links for plaquette action density (matches openQCD) --- src/YM/YMflow.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 42ed545..2ba3d23 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -303,6 +303,10 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T,G,NN}, lp::SpaceParm{N,M,B,D}, ymws: for it in 1:lp.iL[end] Eslc[it,ipl] = 2*Etmp[it] end + if OBC ## Spatial plaquettes at time boundary count half (Luescher) + Eslc[1,ipl] = Etmp[1] + Eslc[end,ipl] = Etmp[end] + end end end @@ -403,6 +407,12 @@ function Eoft_clover(Eslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMwor for it in 1:lp.iL[end] Eslc[it,ipl1] = Etmp[it]/8 end + + ## tentative (weight 1/2 on spatial planes on boundary) + #if (B == BC_OPEN) && (ipl1 >= 4) + # Eslc[end,ipl1] = Etmp[end]/16 + # Eslc[1,ipl1] = Etmp[1]/16 + #end CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_et!(ymws.rm, ymws.frc2, lp) @@ -412,6 +422,12 @@ function Eoft_clover(Eslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMwor Eslc[it,ipl2] = Etmp[it]/8 end + ## tentative (weight 1/2 on spatial planes on boundary) + #if (B == BC_OPEN) && (ipl2 >= 4) + # Eslc[end,ipl2] = Etmp[end]/16 + # Eslc[1,ipl2] = Etmp[1]/16 + #end + return nothing end From 58685f63a243c06ab0a3b7209ccbc3225b4c448c Mon Sep 17 00:00:00 2001 From: Nicolas Lang Date: Tue, 20 May 2025 16:58:19 +0200 Subject: [PATCH 4/5] changed time-boundary of Eoft_clover and Eoft_plaq to match openQCD --- src/YM/YMflow.jl | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 8004653..92c8e2c 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -430,29 +430,25 @@ function Eoft_clover(Eslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMwor end Etmp .= reshape(Array(CUDA.reduce(+, ymws.rm;dims=tp)),lp.iL[end])/V3 for it in 1:lp.iL[end] - Eslc[it,ipl1] = Etmp[it]/8 + if (B == BC_OPEN) && (it == 1 || it == lp.iL[end]) + Eslc[it,ipl1] = 0.0 + else + Eslc[it,ipl1] = Etmp[it]/8 + end end - ## tentative (weight 1/2 on spatial planes on boundary) - #if (B == BC_OPEN) && (ipl1 >= 4) - # Eslc[end,ipl1] = Etmp[end]/16 - # Eslc[1,ipl1] = Etmp[1]/16 - #end - CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_et!(ymws.rm, ymws.frc2, lp) end Etmp .= reshape(Array(CUDA.reduce(+, ymws.rm;dims=tp)),lp.iL[end])/V3 for it in 1:lp.iL[end] - Eslc[it,ipl2] = Etmp[it]/8 + if (B == BC_OPEN) && (it == 1 || it == lp.iL[end]) + Eslc[it,ipl2] = 0.0 + else + Eslc[it,ipl2] = Etmp[it]/8 + end end - ## tentative (weight 1/2 on spatial planes on boundary) - #if (B == BC_OPEN) && (ipl2 >= 4) - # Eslc[end,ipl2] = Etmp[end]/16 - # Eslc[1,ipl2] = Etmp[1]/16 - #end - return nothing end From 9ba07b506aa4d8919485692be266c81bb5eefd36 Mon Sep 17 00:00:00 2001 From: Nicolas Lang Date: Mon, 26 May 2025 12:48:04 +0200 Subject: [PATCH 5/5] fixed the OBC branch in the Zeuthen term - now results match openQCD --- src/YM/YM.jl | 2 +- src/YM/YMflow.jl | 30 ++++++++++++++++++++++++++---- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 5578730..8f862a4 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -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! diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 92c8e2c..62abde1 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -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