Skip to content

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

Added DiracParam struct and changed functions for its inclussion

parent 14c116f7
Hide whitespace changes
Inline Side-by-side
Showing with 272 additions and 16 deletions
src/Dirac/Dirac.jl
View file @ 81c60e34
......@@ -19,6 +19,15 @@ using ..Fields
using ..YM
using ..Spinors
struct DiracParam{T}
rep
m0::T
csw::T
th::NTuple{4,T}
rw::T
end
struct DiracWorkspace{T}
sr
sp
......@@ -34,27 +43,27 @@ struct DiracWorkspace{T}
return new{T}(sr,sp,sAp,st)
end
end
export DiracWorkspace
export DiracWorkspace, DiracParam
function Dw!(so, U, si, m0, lp::SpaceParm)
function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm)
@timeit "Dw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, m0, [1,1,1,1], lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th , lp)
end
end
return nothing
end
function DwdagDw!(so, U, si, m0, st, lp::SpaceParm)
function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm)
@timeit "DwdagDw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(st, U, si, m0, [1,1,1,1], lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(st, U, si, dpar.m0, dpar.th, lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, st, m0, [1,1,1,1], lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, st, dpar.m0, dpar.th, lp)
end
end
......@@ -62,6 +71,17 @@ function DwdagDw!(so, U, si, m0, st, lp::SpaceParm)
end
function g5Dw!(so, U, si, dpar, lp::SpaceParm)
@timeit "DwdagDw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, lp)
end
end
return nothing
end
function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
b = Int64(CUDA.threadIdx().x); r = Int64(CUDA.blockIdx().x)
......@@ -71,13 +91,22 @@ function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
# - Spinor can be periodic as long as 0 at x_0=0
@inbounds begin
so[b,r] = (4+m0)*si[b,r]
for id in 1:4
bu, ru = up((b,r), id, lp)
bd, rd = dw((b,r), id, lp)
bu1, ru1 = up((b,r), 1, lp)
bd1, rd1 = dw((b,r), 1, lp)
bu2, ru2 = up((b,r), 2, lp)
bd2, rd2 = dw((b,r), 2, lp)
bu3, ru3 = up((b,r), 3, lp)
bd3, rd3 = dw((b,r), 3, lp)
bu4, ru4 = up((b,r), 4, lp)
bd4, rd4 = dw((b,r), 4, lp)
so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
end
so[b,r] -= 0.5*( th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) +
th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) +
th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) +
th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) )
end
return nothing
......@@ -89,19 +118,246 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D}
@inbounds begin
so[b,r] = (4+m0)*si[b,r]
bu1, ru1 = up((b,r), 1, lp)
bd1, rd1 = dw((b,r), 1, lp)
bu2, ru2 = up((b,r), 2, lp)
bd2, rd2 = dw((b,r), 2, lp)
bu3, ru3 = up((b,r), 3, lp)
bd3, rd3 = dw((b,r), 3, lp)
bu4, ru4 = up((b,r), 4, lp)
bd4, rd4 = dw((b,r), 4, lp)
so[b,r] -= 0.5*( th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) +
th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) +
th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) +
th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) )
so[b,r] = dmul(Gamma{5}, so[b,r])
end
return nothing
end
############################### HMC for fermions ###################################
function pfrandomize!(f,lp::SpaceParm,dpar::DiracParam,t::Int64=0) #DiracParam to apply dirac operator later?
if dpar.rep == SU3fund && lp.ndim == 4
@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
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_pf!(f,p,lp,t)
end
end
end
return nothing
end
function krnl_assign_pf!(f::AbstractArray{T}, p ,lp::SpaceParm, t::Int64) where {T} #only valid for SU3fund for now. Check performance and maybe change it for the tuple gen in the krnl
@inbounds begin
b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x)
if t == 0
f[b,r] = Spinor(map(x->SU3fund(x[b,1,r,1] + im* x[b,1,r,2],
x[b,2,r,1] + im* x[b,2,r,2],
x[b,3,r,1] + im* x[b,3,r,2]),p))
elseif point_time((b,r),lp) == t
f[b,r] = Spinor(map(x->SU3fund(x[b,1,r,1] + im* x[b,1,r,2],
x[b,2,r,1] + im* x[b,2,r,2],
x[b,3,r,1] + im* x[b,3,r,2]),p))
end
end
return nothing
end
#=
function fermion_action(U, ymws::YMworkspace, lp::SpaceParm, dws::DiracWorkspace,dpar::DiracParam)
if abs(dws.csw) < 1.0E-10
@timeit "Wilson fermions action" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_pfermion_action!(ymws.cm, dws.sp, U, lp,dpar) #which s??
end
end
else
## Csw improved action add it ?
end
S = CUDA.mapreduce(real, +, ymws.cm)
end
function krnl_pfermion_action!(sad, phi ,U::AbstractArray{T}, lp::SpaceParm, dpar::DiracParam)
b = Int64(CUDA.threadIdx().x);
r = Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
m0 = dpar.m0[1] # only for degenerate pairs of fermions
rw = dpar.rw
th = dpar.th
@inbounds begin
sad[I]] = (4+m0)*Diracdot(phi[b,r],phi[b,r])
for id in 1:4
bu, ru = up((b,r), id, lp)
bd, rd = dw((b,r), id, lp)
so[b,r] -= ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],si[bu,ru]) +
conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],si[bd,rd]) )/2
sad[I] -= Diracdot( phi[b,r] , dmul( Gamma{4}, ( th[id]*gpmul(Pgamma{id,-1},U[b,id,r],phi[bu,ru]) +
conj(th[id])*gdagpmul(Pgamma{id,+1},U[bd,id,rd],phi[bd,rd]) )/2))
end
so[b,r] = dmul(Gamma{5}, so[b,r])
end # Do it for the (DDdag )-1
return nothing
end
function hamiltonian_pf(mom, U, lp, gp, ymws, dws, dpar)
@timeit "Computing Hamiltonian" begin
K = CUDA.mapreduce(norm2, +, mom)/2
V = gauge_action(U, lp, gp, ymws)
F = fermion_action(U, ymws, lp, dws, dpar) # check which dws will have the actual configuration
#Calculo de phi* (DwDwdag)^-1 phi aprovechando (DwDwdag)^-1)phi del calculo para la fuerza,
#hecho antes y guardado en dws. Sumarlo al output. Save te fermion action density of sum it directly?
end
return K+V
end
# 8.44 GL. We need DDwdag not Dwdag Dw ?????
function force_fermion(ffrc, U, dws::DiracWorkspace, lp::SpaceParm)
if abs(dws.csw) < 1.0E-10
@timeit "Wilson fermions action" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_force_fermion!()
end
end
else
## Csw improved action
end
end
function HMCF!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}, dws::DiracWorkspace{T}; noacc=false) where T
@timeit "HMC trayectory" begin
ymws.U1 .= U
randomize!(ymws.mom, lp, ymws)
# randomizar campo xi y guardarphi=Dw xi
hini = hamiltonian_pf(ymws.mom, U, lp, gp, ymws, dws, dpar)
#modificar la funcion H para que incluya el t'ermino de pseudofermiones (que depende de los U)
MD!(ymws.mom, U, int, lp, gp, ymws)
#aplicar MD con los terminos de fuerza de pseudofermiones
dh = hamiltonian_pf(ymws.mom, U, lp, gp, ymws) - hini
pacc = exp(-dh)
acc = true
if (noacc)
return dh, acc
end
if (pacc < 1.0)
r = rand()
if (pacc < r)
U .= ymws.U1
acc = false
end
end
#probablemente resetar los campos del dws?
U .= unitarize.(U)
end
return dh, acc
end
#modificar para que aplique ambas fuerzas
function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat}
@timeit "MD evolution" begin
ee = int.eps*gp.beta/gp.ng
force_gauge(ymws, U, gp.c0, gp, lp)
mom .= mom .+ (int.r[1]*ee) .* ymws.frc1
for i in 1:int.ns
k = 2
off = 1
for j in 1:NI-1
U .= expm.(U, mom, int.eps*int.r[k])
if k == NI
off = -1
end
k += off
force_gauge(ymws, U, gp.c0, gp, lp)
if (i < int.ns) && (k == 1)
mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1
else
mom .= mom .+ (int.r[k]*ee) .* ymws.frc1
end
if k == NI
off = -1
end
k += off
end
end
end
return nothing
end
export Dw!, DwdagDw!
=#
export Dw!, DwdagDw!, g5Dw!, pfrandomize!
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