From e18183d5177105a171c5d87646fa6c544c94123d Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Sat, 25 Sep 2021 09:35:56 +0200 Subject: [PATCH] Some cleanup --- src/Groups/GroupSU2.jl | 8 +++---- src/YM/YM.jl | 22 +++++++++++++++++++ src/YM/YMact.jl | 48 ++++++++++++++++++++++++------------------ src/YM/YMfields.jl | 2 -- src/YM/YMhmc.jl | 47 ++++++----------------------------------- 5 files changed, 61 insertions(+), 66 deletions(-) diff --git a/src/Groups/GroupSU2.jl b/src/Groups/GroupSU2.jl index 32d85e7..737c872 100644 --- a/src/Groups/GroupSU2.jl +++ b/src/Groups/GroupSU2.jl @@ -31,9 +31,9 @@ Base.one(::Type{SU2{T}}) where T <: AbstractFloat = SU2{T}(one(T),zero(T)) Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2{T}}) where T <: AbstractFloat = exp(SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T))) """ - function normalize(a::SU2) + function normalize(a::T) where {T <: Group} -Return a normalized element of `SU(2)` +Return a normalized element of the group. """ function normalize(a::SU2{T}) where T <: AbstractFloat dr = sqrt(abs2(a.t1) + abs2(a.t2)) @@ -80,7 +80,7 @@ function isgroup(a::SU2{T}) where T <: AbstractFloat end """ - function Base.exp(a::SU2alg, t::Number=1) + function Base.exp(a::T, t::Number=1) where {T <: Algebra} Computes `exp(a)` """ @@ -120,7 +120,7 @@ end """ - function expm(g::SU2, a::SU2alg) + function expm(g::G, a::A) where {G <: Algebra, A <: Algebra} Computes `exp(a)*g` diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 6ae1f27..8bcb7d5 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -16,12 +16,21 @@ using CUDA, Random, StructArrays using ..Space using ..Groups +import Base.show + struct GaugeParm{T} beta::T cG::NTuple{2,T} ng::Int64 end export GaugeParm +function Base.show(io::IO, gp::GaugeParm) + + println(io, "beta: ", gp.beta) + println(io, "Ngauge: ", gp.beta) + + return nothing +end struct YMworkspace{T} GRP @@ -70,6 +79,19 @@ struct YMworkspace{T} end end export YMworkspace +function Base.show(io::IO, ymws::YMworkspace) + + println(io, "Workspace for Group: ", ymws.GRP) + println(io, " Algebra: ", ymws.ALG) + println(io, "Precision: ", ymws.PRC) + if ymws.fpln == nothing + println(io, " - Running in memory efficient mode") + else + println(io, " - Running in computing efficient mode") + end + return nothing +end + include("YMfields.jl") export field, field_pln, randomize!, zero!, norm2 diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index 8302777..6e07026 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -9,29 +9,31 @@ ### created: Mon Jul 12 18:31:19 2021 ### -function krnl_plaq!(plx, U, ipl, lp::SpaceParm) - - id1, id2 = lp.plidx[ipl] - b, r = CUDA.threadIdx().x, CUDA.blockIdx().x - bu1, ru1 = up((b, r), id1, lp) - bu2, ru2 = up((b, r), id2, lp) - - @inbounds plx[b, r] = tr(U[b,id1,r]*U[bu1,id2,ru1] / (U[b,id2,r]*U[bu2,id1,ru2])) - - return nothing -end - -function krnl_plaq!(plx, U, lp::SpaceParm) +function krnl_plaq!(plx, U, lp::SpaceParm{N,M,D}) where {T,N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x + + Ush = @cuStaticSharedMem(T, (D,2)) + plx[b,r] = zero(plx[b,r]) - @inbounds for ipl in 1:lp.npls - id1, id2 = lp.plidx[ipl] - - bu1, ru1 = up((b, r), id1, lp) - bu2, ru2 = up((b, r), id2, lp) - - plx[b,r] += tr(U[b,id1,r]*U[bu1,id2,ru1] / (U[b,id2,r]*U[bu2,id1,ru2])) + for id1 in 1:N-1 + if ru2 == r + gt2 = Ush[bu2,1] + else + gt2 = U[bu2,id1,ru2] + end + for id2 = id1+1:N + if ru1 == r + gt1 = Ush[bu1,2] + else + gt1 = U[bu1,id2,ru1] + end + sync_threads() + + bu1, ru1 = up((b, r), id1, lp) + bu2, ru2 = up((b, r), id2, lp) + + plx[b,r] += tr(Ush[b,1]*gt1 / (Ush[b,2]*gt2)) end return nothing @@ -139,6 +141,12 @@ function krnl_add_force_plns!(frc::AbstractArray{T}, fpl, lp::SpaceParm{N,M,D}) return nothing end +""" + function force_wilson(ymws::YMworkspace, U, lp::SpaceParm) + +Computes the force deriving from the Wilson plaquette action, without +the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` +""" function force_wilson(ymws::YMworkspace, U, lp::SpaceParm) if ymws.fpln == nothing diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl index 6f5cbfa..3797766 100644 --- a/src/YM/YMfields.jl +++ b/src/YM/YMfields.jl @@ -9,8 +9,6 @@ ### created: Thu Jul 15 15:16:47 2021 ### -un(t) = t <: Union{Group, Complex} - function field(::Type{T}, lp::SpaceParm) where {T} sz = lp.bsz, lp.ndim, lp.rsz diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 5301fa0..56fcc68 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -9,7 +9,13 @@ ### created: Thu Jul 15 11:27:28 2021 ### -function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat +""" + + function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) + +Returns the value of the Wilson plaquette action for the configuration U. +""" +function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: A CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, lp) @@ -64,24 +70,6 @@ function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc return dh, acc end -function krnl_updt!(mom::AbstractArray{TF}, frc, eps1, U::AbstractArray{TU}, eps2, lp::SpaceParm{N,M,D}) where {TU,TF, N,M,D} - - b, r = CUDA.threadIdx().x, CUDA.blockIdx().x - - Ush = @cuStaticSharedMem(TU, D) - Fsh = @cuStaticSharedMem(TF, D) - - @inbounds for id in 1:lp.ndim - Ush[b] = U[b,id,r] - Fsh[b] = frc[b,id,r] - - mom[b,id,r] = mom[b,id,r] + eps1 * Fsh[b] - U[b,id,r] = expm(Ush[b], mom[b,id,r], eps2) - end - - return nothing -end - function OMF4!(mom, U, eps, ns, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat r1::T = 0.08398315262876693 @@ -94,46 +82,25 @@ function OMF4!(mom, U, eps, ns, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspa ee = eps*gp.beta/gp.ng force_wilson(ymws, U, lp) for i in 1:ns - # STEP 1 -# CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_updt!(ymws.mom, ymws.frc1, r1*ee, U, eps*r2, lp) -# end mom .= mom .+ (r1*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r2) - # STEP 2 force_wilson(ymws, U, lp) -# CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_updt!(ymws.mom, ymws.frc1, r3*ee, U, eps*r4, lp) -# end mom .= mom .+ (r3*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r4) - # STEP 3 force_wilson(ymws, U, lp) -# CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_updt!(ymws.mom, ymws.frc1, r5*ee, U, eps*r6, lp) -# end mom .= mom .+ (r5*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r6) - # STEP 4 force_wilson(ymws, U, lp) -# CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_updt!(ymws.mom, ymws.frc1, r5*ee, U, eps*r4, lp) -# end mom .= mom .+ (r5*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r4) - # STEP 5 force_wilson(ymws, U, lp) -# CUDA.@sync begin -# CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_updt!(ymws.mom, ymws.frc1, r3*ee, U, eps*r2, lp) -# end mom .= mom .+ (r3*ee) .* ymws.frc1 U .= expm.(U, mom, eps*r2) - # STEP 6 force_wilson(ymws, U, lp) mom .= mom .+ (r1*ee) .* ymws.frc1 end