diff --git a/src/Groups/FundamentalSU2.jl b/src/Groups/FundamentalSU2.jl index 74de503..b0bc666 100644 --- a/src/Groups/FundamentalSU2.jl +++ b/src/Groups/FundamentalSU2.jl @@ -10,15 +10,18 @@ ### SU2fund(a::T, b::T) where T <: AbstractFloat = SU2fund{T}(complex(a), complex(b)) -dag(a::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(conj(b.t1), -b.t2) +dag(a::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(conj(a.t1), -a.t2) norm(a::SU2fund{T}) where T <: AbstractFloat = sqrt((abs2(a.t1) + abs2(a.t2))/2) norm2(a::SU2fund{T}) where T <: AbstractFloat = (abs2(a.t1) + abs2(a.t2))/2 tr(g::SU2fund{T}) where T <: AbstractFloat = complex(real(g.t1), 0.0) -dot(g1::SU2fund{T},g2::SU2fund{T}) where T <: AbstractFloat = real(conj(g1.t1)*g2.t1+g1.t2*conj(g2.t2)) +dot(g1::SU2fund{T},g2::SU2fund{T}) where T <: AbstractFloat = real(conj(g1.t1)*g2.t1+g1.t2*conj(g2.t2))/2 +projalg(g::SU2fund{T}) where T <: AbstractFloat = SU2alg{T}(imag(g.t2)/2, real(g.t2)/2, imag(g.t1)/2) -Base.:*(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1*b.t1-a.t2*conj(b.t2),a.t1*b.t2+a.t2*conj(b.t1)) -Base.:/(a::SU2fund{T},b::SU2{T}) where T <: AbstractFloat = SU2fund{T}(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1) -Base.:\(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1)) +Base.:*(a::SU2fund{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}((a.t1*b.t1-a.t2*conj(b.t2))/2,(a.t1*b.t2+a.t2*conj(b.t1))/2) +Base.:*(a::SU2fund{T},b::SU2{T}) where T <: AbstractFloat = SU2fund{T}( a.t1*b.t1-a.t2*conj(b.t2) , a.t1*b.t2+a.t2*conj(b.t1)) +Base.:*(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}( a.t1*b.t1-a.t2*conj(b.t2) , a.t1*b.t2+a.t2*conj(b.t1)) +Base.:/(a::SU2fund{T},b::SU2{T}) where T <: AbstractFloat = SU2fund{T}(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1) +Base.:\(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1)) Base.:+(a::SU2fund{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1+b.t1,a.t2+b.t2) Base.:-(a::SU2fund{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1-b.t1,a.t2-b.t2) diff --git a/src/Groups/SU2Types.jl b/src/Groups/SU2Types.jl index 391660c..3e43124 100644 --- a/src/Groups/SU2Types.jl +++ b/src/Groups/SU2Types.jl @@ -37,6 +37,8 @@ Base.zero(::Type{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(zero(T),zero( Base.zero(::Type{M2x2{T}}) where T <: AbstractFloat = M2x2{T}(zero(T),zero(T),zero(T),zero(T)) Base.one(::Type{SU2{T}}) where T <: AbstractFloat = SU2{T}(one(T),zero(T)) Base.one(::Type{M2x2{T}}) where T <: AbstractFloat = M2x2{T}(one(T),zero(T),zero(T),one(T)) +Base.one(::Type{SU2fund{T}}) where T <: AbstractFloat = SU2fund{T}(2*one(T),zero(T)) +Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2fund{T}}) where T <: AbstractFloat = SU2fund{T}(complex(randn(rng,T),randn(rng,T)),complex(randn(rng,T),randn(rng,T))) Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T)) Random.rand(rng::AbstractRNG, ::Random.SamplerType{SU2{T}}) where T <: AbstractFloat = exp(SU2alg{T}(randn(rng,T),randn(rng,T),randn(rng,T))) diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 9c126a5..cae70bf 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -43,7 +43,8 @@ export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3 include("Scalar/Scalar.jl") using .Scalar -export ScalarParam, ScalarWorkspace +export ScalarParm, ScalarWorkspace export scalar_action, force_scalar +export HMC_Scalar!, randomize! end # module diff --git a/src/Scalar/Scalar.jl b/src/Scalar/Scalar.jl index 85dfb0a..c175457 100644 --- a/src/Scalar/Scalar.jl +++ b/src/Scalar/Scalar.jl @@ -14,36 +14,42 @@ module Scalar using CUDA, Random using ..Space using ..Groups +using ..Fields +using ..MD using ..YM +import ..YM: HMC!, randomize!, MD!, force_gauge import Base.show -struct ScalarParam{N,T} +struct ScalarParm{N,T} kap::NTuple{N,T} eta::NTuple{N,T} end -function Base.show(io::IO, sp::ScalarParam{N,T}) where {N,T} +function Base.show(io::IO, sp::ScalarParm{N,T}) where {N,T} println(io, "Number of scalar fields: ", N) print(io, "Kappas: ") for i in 1:N print(io, " ", sp.kap[i]) end - println("\netas: ") + print("\netas: ") for i in 1:N - print(io, " ", sp.kap[i]) + print(io, " ", sp.eta[i]) end println("\n") end +export ScalarParm struct ScalarWorkspace{T} frc1 mom - function ScalarWorkspace(::Type{T}, lp::SpaceParm) where {T <: AbstractFloat} + Phi + function ScalarWorkspace(::Type{T}, n, lp::SpaceParm) where {T <: AbstractFloat} - return ScalarWorkspace(vector_field(SU2fund{T}, lp), - vector_field(SU2fund{T}, lp)) + return new{T}(nscalar_field(SU2fund{T}, n, lp), + nscalar_field(SU2fund{T}, n, lp), + nscalar_field(SU2fund{T}, n, lp)) end end export ScalarWorkspace @@ -54,4 +60,10 @@ export scalar_action include("ScalarForce.jl") export force_scalar +include("ScalarFields.jl") +export randomize! + +include("ScalarHMC.jl") +export HMC! + end diff --git a/src/Scalar/ScalarAction.jl b/src/Scalar/ScalarAction.jl index eeb439f..cb01d88 100644 --- a/src/Scalar/ScalarAction.jl +++ b/src/Scalar/ScalarAction.jl @@ -9,7 +9,7 @@ ### created: Tue Oct 5 11:53:49 2021 ### -function scalar_action(U, Phi, lp::SpaceParm, sp::ScalarParam, ymws::YMworkspace{T}) where {T <: AbstractFloat} +function scalar_action(U, Phi, lp::SpaceParm, sp::ScalarParm, ymws::YMworkspace{T}) where {T <: AbstractFloat} CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_act!(ymws.rm, U, Phi, sp, lp) @@ -19,7 +19,7 @@ function scalar_action(U, Phi, lp::SpaceParm, sp::ScalarParam, ymws::YMworkspace return S end -function krnl_act!(act, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParam{NP,T}, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D} +function krnl_act!(act, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParm{NP,T}, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x @@ -39,14 +39,14 @@ function krnl_act!(act, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::Scalar for id in 1:N bu, ru = up((b, r), id, lp) - if ru == r - Pup = ntuple(i -> Psh[bu,i], NP) - else - Pup = ntuple(i -> Phi[bu,i,ru], NP) - end - for i in 1:NP - act[b,r] += -2*sp.kap[i]*dot(Psh[b,i],Ush[b,id]*Pup[i]) + if ru == r + Pup = Psh[bu,i] + else + Pup = Phi[bu,i,ru] + end + + act[b,r] += -2*sp.kap[i]*dot(Psh[b,i],Ush[b,id]*Pup) end end diff --git a/src/Scalar/ScalarForce.jl b/src/Scalar/ScalarForce.jl index fb3d703..b36b8ad 100644 --- a/src/Scalar/ScalarForce.jl +++ b/src/Scalar/ScalarForce.jl @@ -9,17 +9,17 @@ ### created: Wed Oct 6 15:39:07 2021 ### -function force_scalar(ymws::YMworkspace, sws::ScalarWorkspace, U, Phi, sp::ScalarParam, lp::SpaceParm) +function force_scalar(ymws::YMworkspace, sws::ScalarWorkspace, U, Phi, sp::ScalarParm, gp::GaugeParm, lp::SpaceParm) CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_scalar!(ymws.frc1,sws.frc,U,Phi,sp,lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_scalar!(ymws.frc1,sws.frc1,U,Phi,sp,gp,lp) end return nothing end -function krnl_force_scalar!(fgauge, fscalar, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParam{NP,T}, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D} +function krnl_force_scalar!(fgauge, fscalar, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParm{NP,T}, gp::GaugeParm, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D} b, r = CUDA.threadIdx().x, CUDA.blockIdx().x @@ -28,33 +28,34 @@ function krnl_force_scalar!(fgauge, fscalar, U::AbstractArray{TG}, Phi::Abstract for id in 1:N Ush[b,id] = U[b,id,r] + fgauge[b,id,r] = (gp.beta/gp.ng)*fgauge[b,id,r] end for i in 1:NP Psh[b,i] = Phi[b,i,r] end sync_threads() - for n in 1:NP - fscalar[b,n,u] = zero(TS) + for i in 1:NP + fscalar[b,i,r] = zero(TS) for id in 1:N bu, ru = up((b,r), id, lp) bd, rd = dw((b,r), id, lp) if ru == r - p1 = Ush[b,id]*Psh[bu,n] + p1 = Ush[b,id]*Psh[bu,i] else - p1 = Ush[b,id]*Phi[bu,n,ru] + p1 = Ush[b,id]*Phi[bu,i,ru] end if rd == r - fscalar[b,n,u] += (-2*sp.kap[n])*(p1 + dag(Psh[bd,n])*Ush[bd,id]) + fscalar[b,i,r] += (2*sp.kap[i])*(p1 + dag(Ush[bd,id])*Psh[bd,i]) else - fscalar[b,n,u] += (-2*sp.kap[n])*(p1 + dag(Phi[bd,n,rd])*U[bd,id,rd]) - end + fscalar[b,i,r] += (2*sp.kap[i])*(p1 + dag(U[bd,id,rd])*Phi[bd,i,rd]) + end - fgauge[b,id,r] += (2*sp.kap[n])*projalg(p1*dag(Psh[b,n])) + fgauge[b,id,r] -= (2*sp.kap[i])*projalg(p1*dag(Psh[b,i])) end - fscalar[b,n,r] += ( 2 + 2*sp.eta[n]*(dot(Psh[b,n],Psh[b,n])-1) ) * Psh[b,n] + fscalar[b,i,r] -= ( 2 + 4*sp.eta[i]*(dot(Psh[b,i],Psh[b,i])-1) ) * Psh[b,i] end diff --git a/tests/test_SU2.jl b/tests/test_SU2.jl index 12f3739..ed12143 100644 --- a/tests/test_SU2.jl +++ b/tests/test_SU2.jl @@ -50,3 +50,142 @@ println(mp) println(projalg(mp)) println(projalg(ga)) +println("## HERE test SU(2) fundamental") +a = rand(SU2fund{T}) +println(a) + +ft(x) = LatticeGPU.dot(x,x) +f2(x) = LatticeGPU.tr(dag(x)*x) +f3(x) = LatticeGPU.norm2(x) + +S = ft(a) +println("Check 1: ", S) +S = f2(a) +println("Check 2: ", S) +S = f3(a) +println("Check 3: ", S) + +h = 0.000001 +xh = SU2fund{T}(a.t1+h, a.t2) +Sp = ft(xh) +xh = SU2fund{T}(a.t1-h, a.t2) +Sm = ft(xh) +println("Numerical derivative: ", (Sp-Sm)/(2*h)) +xh = SU2fund{T}(a.t1, a.t2+h) +Sp = ft(xh) +xh = SU2fund{T}(a.t1, a.t2-h) +Sm = ft(xh) +println("Numerical derivative: ", (Sp-Sm)/(2*h)) +h = 0.000001im +xh = SU2fund{T}(a.t1+h, a.t2) +Sp = ft(xh) +xh = SU2fund{T}(a.t1-h, a.t2) +Sm = ft(xh) +println("Numerical derivative: ", 1im * (Sp-Sm)/(2*h)) +xh = SU2fund{T}(a.t1, a.t2+h) +Sp = ft(xh) +xh = SU2fund{T}(a.t1, a.t2-h) +Sm = ft(xh) +println("Numerical derivative: ", 1im * (Sp-Sm)/(2*h)) + +ad = a +println("Exact derivative: ", ad) + + + +println("## Hamiltonian dynamics") +println(" # Mass terms") +H(p,a) = LatticeGPU.norm2(p)/2 + LatticeGPU.norm2(a) +function MD!(p, a, ns, ee) + for i in 1:ns + p = p - ee*a + a = a + ee*p + p = p - ee*a + end + return p, a +end + +a = rand(SU2fund{T}) +p = rand(SU2fund{T}) +println(H(p,a)) + +ee = 0.0001 +ns = 10000 +p, a = MD!(p, a, ns, ee) + +println(H(p,a)) + +println(" # Interaction terms") +H(p,a, M,kap) = LatticeGPU.norm2(p)/2 - (2*kap)*LatticeGPU.dot(a, M) +function MD!(p, a, M, kap, ns, ee) + for i in 1:ns + p = p + 0.5*( (2*kap)*ee*M ) + a = a + ee*p + p = p + 0.5*( (2*kap)*ee*M ) + end + return p, a +end + +a = rand(SU2fund{T}) +p = rand(SU2fund{T}) +M = rand(SU2fund{T}) +kap = 0.5 +println(H(p,a, M,kap)) + +ee = 0.0001 +ns = 10000 +p, a = MD!(p, a, M,kap, ns, ee) + +println(H(p,a, M,kap)) + +println(" # Gauge Interaction terms") +H(p,a, M,kap) = LatticeGPU.norm2(p)/2 - (2*kap)*LatticeGPU.tr(a*M) +function MD!(p, a, M, kap, ns, ee) + for i in 1:ns + p = p + 0.5*( - (2*kap)*ee*projalg(a*M) ) + a = expm(a, p, ee) + p = p + 0.5*( - (2*kap)*ee*projalg(a*M) ) + end + return p, a +end + +a = rand(SU2{T}) +p = rand(SU2alg{T}) +M = rand(SU2fund{T}) +#M = SU2fund{T}(0.0+rand()im, rand()+rand()im) +kap = 1.1 +println(H(p,a, M,kap)) + +ee = 0.0001 +ns = 10000 +p, a = MD!(p, a, M,kap, ns, ee) + +println(H(p,a, M,kap)) + +println(" # Gauge Interaction terms") +H(p,a, M1,M2,kap) = LatticeGPU.norm2(p)/2 - (2*kap)*LatticeGPU.dot(M1,a*M2) +#H(p,a, M1,M2,kap) = LatticeGPU.norm2(p)/2 - (2*kap)*LatticeGPU.tr(a*(M2*dag(M1))) +function MD!(p, a, M1, M2, kap, ns, ee) + for i in 1:ns + p = p + 0.5*( - (2*kap)*ee*projalg(a*M2*dag(M1)) ) + a = expm(a, p, ee) + p = p + 0.5*( - (2*kap)*ee*projalg(a*M2*dag(M1)) ) + end + return p, a +end + +a = rand(SU2{T}) +p = rand(SU2alg{T}) +M1 = rand(SU2fund{T}) +M2 = rand(SU2fund{T}) +#M = SU2fund{T}(0.0+rand()im, rand()+rand()im) +kap = 1.1 +println(H(p,a, M1,M2,kap)) + +ee = 0.0001 +ns = 10000 +p, a = MD!(p, a, M1,M2,kap, ns, ee) + +println(H(p,a, M1,M2,kap)) + +