new observables, ScalarObs.jl

src/Scalar/ScalarObs.jl

@@ -0,0 +1,95 @@
+###
+### "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:    YMact.jl
+### created: Mon Jul 12 18:31:19 2021
+###
+
+"""
+    each instance defines fields to store observables in each lattice point
+"""
+
+struct Obs{T}
+    rho2 #ρ^2
+    lphi #L_\phi
+    lalpha #L_\alpha
+    function Obs(::Type{T}, lp::SpaceParm, sp::ScalarParm, gp::GaugeParm) where {T <: AbstractFloat}
+
+        rho2n   = nscalar_field(Complex{T}, length(sp.kap), lp)
+        lphin   = nscalar_field(Complex{T}, length(sp.kap), lp)
+        lalphan = nscalar_field(Complex{T}, length(sp.kap), lp)
+        return new{T}(rho2n, lphin, lalphan)
+    end
+end
+
+"""
+    computes global observables by calling krnl_obs! and summing
+    for all lattice points
+"""
+
+function updt_obs!(obs::Obs{T}, U, Phi, sp::ScalarParm, lp::SpaceParm) where {T <: AbstractFloat}
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_obs!(obs, U, Phi, sp, lp)
+    end
+
+    #summation of global observables
+    rho2 = CUDA.reduce(+, obs.rho2)
+    lphi = CUDA.reduce(+, obs.lphi)
+    lalpha = CUDA.reduce(+, obs.lalpha)
+    return (rho2,lphi,lalpha)
+end
+
+"""
+    CUDA function to compute the observables defined in the Obs struct
+    for each lattice point
+"""
+
+function krnl_obs!(obs::Obs{T}, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParm{NP,T}, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D}
+
+    #thread/block coordinate
+    b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
+
+    Ush = @cuStaticSharedMem(TG, (D,N)) #Links - array (lattice points)x(spacetime dimension)
+    Psh = @cuStaticSharedMem(TS, (D,NP)) #Phield - array (lattice points)x(scalalr dimension)
+
+
+    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()
+
+    #compute obs
+    for i in 1:NP
+        obs.rho2[b,i,r] = norm2( Psh[b,i] )
+        obs.lphi[b,i,r] = zero( obs.lphi[b,i,r] )
+        obs.lalpha[b,i,r] = zero( obs.lalpha[b,i,r] )
+        norm = norm( Psh[b,i] )
+
+        for mu in 1:N
+
+            #up fields
+            bu, ru = up((b, r), id, lp) #thread/block coordinate of up point
+            if (ru == r)
+                phiup = Psh[bu,i]
+            else
+                phiup = Phi[bu,i,ru]
+            end
+
+            normup = norm(phiup)
+
+            obs.lphi[b,i,r] += tr( Psh[b,i], Ush[b,i]*phiup )
+            obs.lalpha[b,i,r] += tr( Psh[b,i]/norm, Ush[b,i]*phiup/normup )
+        end
+    end
+
+    return nothing
+end
+