Added basic support for simulating Scalar field theories

src/Groups/Groups.jl

@@ -25,11 +25,12 @@ export Group, Algebra
 # SU(2) and 2x2 matrix operations
 ##
 include("SU2Types.jl")
-export SU2, SU2alg, M2x2
+export SU2, SU2alg, SU2fund, M2x2
 
 include("GroupSU2.jl")
 include("M2x2.jl")
 include("AlgebraSU2.jl")
+include("FundamentalSU2.jl")
 ## END SU(2)
 
 ##

src/Groups/SU2Types.jl

@@ -27,10 +27,16 @@ struct SU2alg{T} <: Algebra
     t3::T
 end
 
-Base.zero(::Type{SU2alg{T}}) where T <: AbstractFloat = SU2alg{T}(zero(T),zero(T),zero(T))
+struct SU2fund{T}
-Base.zero(::Type{M2x2{T}})   where T <: AbstractFloat = M2x2{T}(zero(T),zero(T),zero(T),zero(T))
+    t1::Complex{T}
-Base.one(::Type{SU2{T}})     where T <: AbstractFloat = SU2{T}(one(T),zero(T))
+    t2::Complex{T}
-Base.one(::Type{M2x2{T}})    where T <: AbstractFloat = M2x2{T}(one(T),zero(T),zero(T),one(T))
+end    
+
+Base.zero(::Type{SU2fund{T}}) where T <: AbstractFloat = SU2fund{T}(zero(T),zero(T))
+Base.zero(::Type{SU2alg{T}})  where T <: AbstractFloat = SU2alg{T}(zero(T),zero(T),zero(T))
+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))
 
 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)))

src/LatticeGPU.jl

@@ -16,7 +16,7 @@ include("Groups/Groups.jl")
 
 using .Groups
 export Group, Algebra
-export SU2, SU2alg, SU3, SU3alg, M3x3, M2x2, U1, U1alg
+export SU2, SU2alg, SU2fund, SU3, SU3alg, M3x3, M2x2, U1, U1alg
 export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat
 
 include("Space/Space.jl")
@@ -33,4 +33,7 @@ export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zer
 export gauge_action, hamiltonian, plaquette, HMC!, OMF4!
 export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3
 
+include("Scalar/Scalar.jl")
+using .Scalar
+
 end # module

src/Scalar/Scalar.jl

@@ -0,0 +1,38 @@
+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    Scalar.jl
+### created: Tue Oct  5 11:53:31 2021
+###                               
+
+module Scalar
+
+using CUDA, Random
+using ..Space
+using ..Groups
+using ..YM
+
+import Base.show
+
+struct ScalarParam{N,T}
+    kap::Ntuple{N,T}
+    eta::Ntuple{N,T}
+end
+function Base.show(io::IO, sp::ScalarParam{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:   ")
+    for i in 1:N
+        print(io, " ", sp.kap[i])
+    end
+    println("\n")
+
+end

src/Scalar/ScalarAction.jl

@@ -0,0 +1,60 @@
+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    ScalarAction.jl
+### created: Tue Oct  5 11:53:49 2021
+###                               
+
+function scalar_action(U, Phi, lp::SpaceParm, sp::ScalarParm, gp::GaugeParm{T}, ymws::YMworkspace{T})
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_act!(ymws.rm, U, Phi, sp, lp)
+    end
+
+    S = CUDA.reduce(+, ymws.rm)
+    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}
+
+
+    b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
+
+    Ush = @cuStaticSharedMem(TG, (D,N))
+    Psh = @cuStaticSharedMem(TS, (D,NP))
+    
+    for id in 1:N
+        Ush[b,id] = U[b,id,r]
+    end
+    for i in 1:NP
+        Psh[b,i] = Phi[b,i,r]
+    end
+    sync_threads()
+    
+    act[b,r] = zero(act[b,r])
+    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])
+        end
+    end
+    
+    for i in 1:NP
+        sdot = dot(Psh[b,i],Psh[b,i])
+        act[b,r] += sdot + sp.eta[i]*(sdot - 1)^2
+    end
+    
+    return nothing
+end
+