mirror of
https://igit.ific.uv.es/alramos/latticegpu.jl.git
synced 2025-05-14 19:23:42 +02:00
Working HMC for one scalar coupled to SU(2)
This commit is contained in:
Alberto Ramos
parent
bd88c213d7
commit
4f9ddb7026
7 changed files with 192 additions and 34 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)))
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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])
|
||||
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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue