Skip to content

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

Test one for SF with fermions

parent bb436c9a
Show whitespace changes
Inline Side-by-side
Showing with 86 additions and 187 deletions
src/Dirac/Dirac.jl
View file @ 1a1f56e3
......@@ -24,7 +24,7 @@ struct DiracParam{T}
m0::T
csw::T
th::NTuple{4,T}
rw::T
ct::T
end
......@@ -45,19 +45,62 @@ struct DiracWorkspace{T}
end
export DiracWorkspace, DiracParam
function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm)
function Dw!(so, U, si, dpar::DiracParam, lp::SpaceParm{4,6,B,D}) where {B,D}
if B == BC_SF_AFWB || B == BC_SF_ORBI
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(si,so,1.0,lp)
end
@timeit "Dw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th , lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, lp)
end
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(so,si,dpar.ct,lp)
end
else
@timeit "Dw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, lp)
end
end
end
return nothing
end
function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm)
function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm{4,6,B,D}) where {B,D}
if B == BC_SF_AFWB || B == BC_SF_ORBI
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(si,so,1.0,lp)
end
@timeit "DwdagDw" begin
CUDA.@sync begin
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_sfbndfix!(st,si,ct,lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, st, dpar.m0, dpar.th, lp)
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(so,st,ct,lp)
end
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(so,si,dpar.ct,lp)
end
else
@timeit "DwdagDw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(st, U, si, dpar.m0, dpar.th, lp)
......@@ -66,6 +109,7 @@ function DwdagDw!(so, U, si, dpar::DiracParam, st, lp::SpaceParm)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, st, dpar.m0, dpar.th, lp)
end
end
end
return nothing
end
......@@ -73,15 +117,44 @@ end
function g5Dw!(so, U, si, dpar, lp::SpaceParm)
@timeit "DwdagDw" begin
if B == BC_SF_AFWB || B == BC_SF_ORBI
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(si,so,1.0,lp)
end
@timeit "Dw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th , lp)
end
end
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(so,si,dpar.ct,lp)
end
else
@timeit "Dw" begin
CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, lp)
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th , lp)
end
end
end
return nothing
end
function krnl_sfbndfix!(so,si,ct,lp::SpaceParm)
b=Int64(CUDA.threadIdx().x)
r=Int64(CUDA.blockIdx().x)
if (point_time((b,r),lp) == 1)
so[b,r] = 0.0*so[b,r]
elseif (point_time((b,r),lp) == 2) || (point_time((b,r),lp) == lp.iL[4])
so[b,r] += (ct-1.0)si[b,r]
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)
......@@ -184,180 +257,6 @@ function krnl_assign_pf!(f::AbstractArray{T}, p ,lp::SpaceParm, t::Int64) where
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)
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
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!, 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