From 3f8ba58848977e81071b475385d8099209b3f270 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Wed, 6 Oct 2021 09:54:24 +0200 Subject: [PATCH] Added basic support for simulating Scalar field theories --- src/Groups/Groups.jl | 3 +- src/Groups/SU2Types.jl | 14 ++++++--- src/LatticeGPU.jl | 5 +++- src/Scalar/Scalar.jl | 38 ++++++++++++++++++++++++ src/Scalar/ScalarAction.jl | 60 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 114 insertions(+), 6 deletions(-) create mode 100644 src/Scalar/Scalar.jl create mode 100644 src/Scalar/ScalarAction.jl diff --git a/src/Groups/Groups.jl b/src/Groups/Groups.jl index b049cd9..bdb3a69 100644 --- a/src/Groups/Groups.jl +++ b/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) ## diff --git a/src/Groups/SU2Types.jl b/src/Groups/SU2Types.jl index d36ec33..391660c 100644 --- a/src/Groups/SU2Types.jl +++ b/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)) -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)) +struct SU2fund{T} + t1::Complex{T} + t2::Complex{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))) diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 5ca9987..cc34475 100644 --- a/src/LatticeGPU.jl +++ b/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 diff --git a/src/Scalar/Scalar.jl b/src/Scalar/Scalar.jl new file mode 100644 index 0000000..5c4b34a --- /dev/null +++ b/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. +### +### 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 diff --git a/src/Scalar/ScalarAction.jl b/src/Scalar/ScalarAction.jl new file mode 100644 index 0000000..c2e80a7 --- /dev/null +++ b/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. +### +### 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 +