From 334ddf813a28888b2c43aea5f67b61f61e1c0794 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Thu, 30 Sep 2021 12:59:00 +0200 Subject: [PATCH] Added routines for Zeuthen flow --- src/Groups/GroupSU3.jl | 46 ++++++++++++++-- src/YM/YMflow.jl | 121 +++++++++++++++++++++++++++-------------- 2 files changed, 121 insertions(+), 46 deletions(-) diff --git a/src/Groups/GroupSU3.jl b/src/Groups/GroupSU3.jl index 266ec1d..d049206 100644 --- a/src/Groups/GroupSU3.jl +++ b/src/Groups/GroupSU3.jl @@ -111,7 +111,7 @@ function projalg(a::SU3{T}) where T <: AbstractFloat sr3ov2::T = 0.866025403784438646763723170752 - ditr = ( imag(a.u11) + imag(a.u22) + 2.0*imag(a.u11*a.u22 - a.u12*a.u21) )/3.0 + ditr = ( imag(a.u11) + imag(a.u22) - 2.0*imag(a.u11*a.u22 - a.u12*a.u21) )/3.0 m12 = (a.u12 - conj(a.u21))/2.0 m13 = (a.u13 - (a.u12*a.u23 - a.u13*a.u22) )/2.0 m23 = (a.u23 - (a.u13*a.u21 - a.u11*a.u23) )/2.0 @@ -122,6 +122,21 @@ function projalg(a::SU3{T}) where T <: AbstractFloat sr3ov2*(ditr)) end +function projalg(a::M3x3{T}) where T <: AbstractFloat + + sr3ov2::T = 0.866025403784438646763723170752 + + ditr = ( imag(a.u11) + imag(a.u22) + 2.0*imag(a.u33) )/3.0 + m12 = (a.u12 - conj(a.u21))/2.0 + m13 = (a.u13 - conj(a.u31))/2.0 + m23 = (a.u23 - conj(a.u32))/2.0 + + return SU3alg{T}(imag( m12 ), imag( m13 ), imag( m23 ), + real( m12 ), real( m13 ), real( m23 ), + (imag(a.u11)-imag(a.u22))/2.0, + sr3ov2*(ditr)) +end + dot(a::SU3alg{T},b::SU3alg{T}) where T <: AbstractFloat = a.t1*b.t1 + a.t2*b.t2 + a.t3*b.t3 + a.t4*b.t4 + a.t5*b.t5 + a.t6*b.t6 + a.t7*b.t7 + a.t8*b.t8 norm2(a::SU3alg{T}) where T <: AbstractFloat = a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2 norm(a::SU3alg{T}) where T <: AbstractFloat = sqrt(a.t1^2 + a.t2^2 + a.t3^2 + a.t4^2 + a.t5^2 + a.t6^2 + a.t7^2 + a.t8^2) @@ -164,9 +179,9 @@ Base.:*(a::M3x3{T},b::M3x3{T}) where T <: AbstractFloat = M3x3{T}(a.u11*b.u11 + function Base.:*(a::SU3{T},b::M3x3{T}) where T <: AbstractFloat - a.u31 = conj(a.u12*a.u23 - a.u13*a.u22) - a.u32 = conj(a.u13*a.u21 - a.u11*a.u23) - a.u33 = conj(a.u11*a.u22 - a.u12*a.u21) + au31 = conj(a.u12*a.u23 - a.u13*a.u22) + au32 = conj(a.u13*a.u21 - a.u11*a.u23) + au33 = conj(a.u11*a.u22 - a.u12*a.u21) return M3x3{T}(a.u11*b.u11 + a.u12*b.u21 + a.u13*b.u31, a.u11*b.u12 + a.u12*b.u22 + a.u13*b.u32, @@ -214,6 +229,24 @@ function Base.:/(a::M3x3{T},b::SU3{T}) where T <: AbstractFloat a.u31*(bu31) + a.u32*(bu32) + a.u33*(bu33)) end +function Base.:\(a::SU3{T},b::M3x3{T}) where T <: AbstractFloat + + au31 = (a.u12*a.u23 - a.u13*a.u22) + au32 = (a.u13*a.u21 - a.u11*a.u23) + au33 = (a.u11*a.u22 - a.u12*a.u21) + + return M3x3{T}(conj(a.u11)*b.u11 + conj(a.u21)*b.u21 + (au31)*b.u31, + conj(a.u11)*b.u12 + conj(a.u21)*b.u22 + (au31)*b.u32, + conj(a.u11)*b.u13 + conj(a.u21)*b.u23 + (au31)*b.u33, + conj(a.u12)*b.u11 + conj(a.u22)*b.u21 + (au32)*b.u31, + conj(a.u12)*b.u12 + conj(a.u22)*b.u22 + (au32)*b.u32, + conj(a.u12)*b.u13 + conj(a.u22)*b.u23 + (au32)*b.u33, + conj(a.u13)*b.u11 + conj(a.u23)*b.u21 + (au33)*b.u31, + conj(a.u13)*b.u12 + conj(a.u23)*b.u22 + (au33)*b.u32, + conj(a.u13)*b.u13 + conj(a.u23)*b.u23 + (au33)*b.u33) + +end + Base.:*(a::Number,b::M3x3{T}) where T <: AbstractFloat = M3x3{T}(a*b.u11, a*b.u12, a*bu13, a*b.u21, a*b.u22, a*bu23, a*b.u31, a*b.u32, a*bu33) @@ -259,6 +292,11 @@ function alg2mat(a::SU3alg{T}) where T <: AbstractFloat return M3x3{T}(u11,u12,u13, u21,u22,u23, u31,u32,u33) end +Base.:*(a::SU3alg,b::SU3) = alg2mat(a)*b +Base.:*(a::SU3,b::SU3alg) = a*alg2mat(b) +Base.:/(a::SU3alg,b::SU3) = alg2mat(a)/b +Base.\:(a::SU3,b::SU3alg) = a\alg2mat(b) + @inline function exp_iter(dch::Complex{T}, tch::T) where T <: AbstractFloat c::NTuple{22, T} = ( 1.957294106339126128e-20, 4.110317623312164853e-19, diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 932f676..c848e7a 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -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)