Skip to content

L
LGPU FPP
Commit 3bd8f080 authored by Fernando Pérez Panadero's avatar Fernando Pérez Panadero
Browse files
Options

Added ct to propagators. Changed src by dws.sr. Added sfcf file

parent f4b15838
Hide whitespace changes
Inline Side-by-side
Showing with 146 additions and 17 deletions
src/Dirac/Dirac.jl
View file @ 3bd8f080
......@@ -51,6 +51,7 @@ struct DiracWorkspace{T}
sAp
st
csw
cs
function DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D}) where {G,T <: AbstractFloat, B,D}
......@@ -60,6 +61,7 @@ struct DiracWorkspace{T}
st = scalar_field(Spinor{4,G}, lp)
csw = tensor_field(U3alg{T},lp)
cs = scalar_field_point(Complex{Float64}, lp)
return new{T}(sr,sp,sAp,st,csw)
end
end
......@@ -436,7 +438,7 @@ end
function pfrandomize!(f,lp::SpaceParm,dpar::DiracParam{T,SU3fund},t::Int64=0) where {T}
function pfrandomize!(f,lp::SpaceParm,t::Int64=0) where {T}
@timeit "Randomize pseudofermion field" begin
p = ntuple(i->CUDA.randn(Float64, lp.bsz, 3, lp.rsz,2),4) # complex generation not suported for Julia 1.5.4
......
src/LatticeGPU.jl
View file @ 3bd8f080
......@@ -64,5 +64,6 @@ include("Solvers/Solvers.jl")
using .Solvers
export CG!
export propagator, bndpropagator
export fP!, fA!, kV!, lV!
end # module
src/Solvers/Propagators.jl
View file @ 3bd8f080
......@@ -18,20 +18,21 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm,
end
pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
src = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
CUDA.@allowscalar src[point_index(CartesianIndex{lp.ndim}(y),lp)...] = Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4))
fill!(dws.sr,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
#src = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
CUDA.@allowscalar dws.sr[point_index(CartesianIndex{lp.ndim}(y),lp)...] = Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sr)
end
g5Dw!(pro,U,src,dpar,lp)
g5Dw!(pro,U,dws.sr,dpar,dws,lp)
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
return pro
end
function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T}
function krnlg5!(src)
......@@ -42,17 +43,16 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm,
end
pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
src = scalar_field(Spinor{4,SU3fund{Float64}},lp)
pfrandomize!(src,lp,dpar,time)
pfrandomize!(dws.sr,lp,dpar,time)
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sr)
end
g5Dw!(pro,U,src,dpar,lp)
g5Dw!(pro,U,dws.sr,dpar,dws,lp)
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sr)
end
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
......@@ -63,8 +63,8 @@ end
function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
Returns the propagator from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'.
Returns the propagator from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not.
For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
"""
function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
......@@ -89,19 +89,64 @@ function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpacePar
end
pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
src = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
fill!(dws.sr,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(src, U, lp, c, s)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(dws.sr, U, lp, c, s)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sr)
end
g5Dw!(pro,U,src,dpar,lp)
g5Dw!(pro,U,dpar.ct*dws.sr,dpar,dws,lp)
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
return pro
end
"""
function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
Returns the propagator from the t=T boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not.
For the propagator from t=00 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64)
"""
function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
function krnlg5!(src)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
src[b,r] = dmul(Gamma{5},src[b,r])
return nothing
end
function krnl_assign_bndsrc!(src,U,lp::SpaceParm, c::Int64, s::Int64)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
if (point_time((b,r),lp) == lp.iL[end])
src[b,r] = gpmul(Pgamma{4,-1},U[b,4,r],Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4)))/2
end
return nothing
end
pro = scalar_field(Spinor{4,SU3fund{Float64}},lp)
fill!(dws.sr,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_bndsrc!(dws.sr, U, lp, c, s)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sr)
end
g5Dw!(pro,U,dpar.ct*dws.sr,dpar,dws,lp)
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol)
return pro
end
src/Solvers/Solvers.jl
View file @ 3bd8f080
......@@ -25,5 +25,8 @@ export CG!
include("Propagators.jl")
export propagator, bndpropagator
include("sfcf.jl")
export fP!, fA!, kV!, lV!
end
src/Solvers/sfcf.jl 0 → 100644
View file @ 3bd8f080
function fP!(fp, dws::DiracWorkspace, prop::AbstractArray,lp::SpaceParm)
function krnl_fP!(sum,pro)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
sum[I] = norm2(pro[b,r])/2
return nothing
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_fP!(dws.cs,prop)
end
fp += reshape(Array(CUDA.reduce(+, dws.cs; dims=(1,2,3))),lp.iL[end])
return nothing
end
function fA!(fa, dws::DiracWorkspace, prop::AbstractArray,lp::SpaceParm)
function krnl_fA!(sum,pro)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
sum[I] = -dot(pro[b,r],dmul(Gamma{4},pro[b,r]))/2
return nothing
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_fA!(dws.cs,prop)
end
fa += reshape(Array(CUDA.reduce(+, dws.cs; dims=(1,2,3))),lp.iL[end])
return nothing
end
function kV!(kv, dws::DiracWorkspace, prop::AbstractArray,lp::SpaceParm)
function krnl_kV!(sum,pro)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
sum[I] = -(norm2(dmul(Gamma{6},pro[b,r])) + norm2(dmul(Gamma{7},pro[b,r])) + norm2(dmul(Gamma{8},pro[b,r])))/6
return nothing
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_kV!(dws.cs,prop)
end
kv += reshape(Array(CUDA.reduce(+, dws.cs; dims=(1,2,3))),lp.iL[end])
return nothing
end
function lV!(lv, dws::DiracWorkspace, prop::AbstractArray,lp::SpaceParm)
function krnl_lV!(sum,pro)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
sum[I] = im*(dot(dmul(Gamma{6},pro[b,r]),dmul(Gamma{5},dmul(Gamma{10},pro[b,r])))
+ dot(dmul(Gamma{7},pro[b,r]),dmul(Gamma{5},dmul(Gamma{11},pro[b,r])))
+ dot(dmul(Gamma{8},pro[b,r]),dmul(Gamma{5},dmul(Gamma{12},pro[b,r]))))/6
return nothing
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_lV!(dws.cs,prop)
end
lv += reshape(Array(CUDA.reduce(+, dws.cs; dims=(1,2,3))),lp.iL[end])
return nothing
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment