From d0c463637fbc629532d367daaac92435b466e890 Mon Sep 17 00:00:00 2001
From: Guilherme telo <gtelo@ific.uv.es>
Date: Tue, 19 Oct 2021 18:23:04 +0200
Subject: [PATCH] new observables, ScalarObs.jl

---
 src/LatticeGPU.jl       |  4 ++
 src/Scalar/ScalarObs.jl | 95 +++++++++++++++++++++++++++++++++++++++++
 src/main/test_scalar.jl | 15 ++++++-
 3 files changed, 113 insertions(+), 1 deletion(-)
 create mode 100644 src/Scalar/ScalarObs.jl

diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index cae70bf..489f78a 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -46,5 +46,9 @@ using .Scalar
 export ScalarParm, ScalarWorkspace
 export scalar_action, force_scalar
 export HMC_Scalar!, randomize!
+#added
+export Obs
+export updt_obs!
+
 
 end # module
diff --git a/src/Scalar/ScalarObs.jl b/src/Scalar/ScalarObs.jl
new file mode 100644
index 0000000..648ab9b
--- /dev/null
+++ b/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
+
diff --git a/src/main/test_scalar.jl b/src/main/test_scalar.jl
index ba98ce3..bcee6f3 100644
--- a/src/main/test_scalar.jl
+++ b/src/main/test_scalar.jl
@@ -10,7 +10,7 @@ lp = SpaceParm{4}((64,64,64,64), (4,4,4,4))
 gp = GaugeParm(6.0, 1.0, (0.0,0.0), 2)
 sp = ScalarParm((0.2,0.3), (1.0,0.4))
 
-NSC = length(sp.kap)
+NSC = length(sp.kap) #number of scalars = # of k coupling
 println("Space  Parameters: ", lp)
 println("Gauge  Parameters: ", gp)
 println("Scalar Parameters: ", sp)
@@ -38,6 +38,10 @@ println("Allocating scalar field")
 Phi = nscalar_field(SCL{PREC}, NSC, lp)
 fill!(Phi, zero(SCL{PREC}))
 
+#initialize observables
+watch = Obs( PREC, lp, sp, gp )
+
+
 println("Initial Action: ")
 @time S = gauge_action(U, lp, gp, ymws) + scalar_action(U, Phi, lp, sp, ymws)
 
@@ -47,6 +51,10 @@ ns  = 20
 
 println("## Thermalization")
 pl = Vector{Float64}()
+rho2_v = Vector{Float64}()
+lphi_v = Vector{Float64}()
+lalpha_v = Vector{Float64}()
+
 for i in 1:10
     @time dh, acc = HMC!(U,Phi, dt,ns,lp, gp, sp, ymws, sws, noacc=true)
     println("# HMC: ", acc, " ", dh)
@@ -60,4 +68,9 @@ for i in 1:10
     println("# HMC: ", acc, " ", dh)
     push!(pl, plaquette(U,lp, gp, ymws))
     println("# Plaquette: ", pl[end], "\n")
+    #measure
+    (rho2, lphi, lalpha) = updt_obs!(watch, U, Phi, sp, lp)
+    push!(rho2_v,rho2)
+    push!(lphi_v,lphi)
+    push!(lalpha_v,lalpha)
 end