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" 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/YMfields.jl b/src/YM/YMfields.jl index 95d6777..2b2db44 100644 --- a/src/YM/YMfields.jl +++ b/src/YM/YMfields.jl @@ -9,16 +9,17 @@ ### 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()) Given an algebra field with natural indexing, this routine sets the components to random Gaussian distributed values. If SF boundary conditions are used, the force at the boundaries is set to zero. """ -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 @@ -28,7 +29,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/YMflow.jl b/src/YM/YMflow.jl index bcb6a35..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 @@ -328,6 +350,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 @@ -426,15 +452,23 @@ 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 - + 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 return nothing diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index d49a8cd..2669cec 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -68,18 +68,19 @@ function hamiltonian(mom, U, lp, gp, ymws) return K+V end + """ - HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc=false) + HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc=false, rng=Random.default_rng(), curng=CUDA.default_rng()) Performs a HMC step (molecular dynamics integration and accept/reject step). The configuration `U` is updated and function returns the energy violation and if the configuration was accepted in a tuple. """ -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) @@ -93,7 +94,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 @@ -105,7 +106,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, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace)