Skip to content

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

Fixed input for propagators and test for Csw

parent 1eb6d458
Hide whitespace changes
Inline Side-by-side
Showing with 39 additions and 31 deletions
src/Dirac/Dirac.jl
View file @ 51a3f002
...@@ -34,6 +34,15 @@ struct DiracParam{T,R} ...@@ -34,6 +34,15 @@ struct DiracParam{T,R}
return new{T,SU3fund}(m0,csw,th,ct) return new{T,SU3fund}(m0,csw,th,ct)
end end
end end
function Base.show(io::IO, dpar::DiracParam{T,R}) where {T,R}``
println(io, "Wilson fermions in the: ", R, " representation")
println(io, " - Bare mass: ", dpar.m0," // Kappa = ",0.5/(dpar.m0+4))
println(io, " - Csw : ", dpar.csw)
println(io, " - c_t: ", dpar.ct)
println(io, " - Theta: ", dpar.th)
return nothing
end
struct DiracWorkspace{T} struct DiracWorkspace{T}
...@@ -353,7 +362,7 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, ztw, lp::SpaceParm{4,M,B ...@@ -353,7 +362,7 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, ztw, lp::SpaceParm{4,M,B
bu2, ru2 = up((b, r), id2, lp) bu2, ru2 = up((b, r), id2, lp)
bd1, rd1 = dw((b, r), id1, lp) bd1, rd1 = dw((b, r), id1, lp)
bd2, rd2 = dw((b, r), id2, lp) bd2, rd2 = dw((b, r), id2, lp)
bdd, rdd = dw((bd1, rd1), id2, lp) bdd, rdd = dw((bd1, rd1), id2, lp)
bud, rud = dw((bu1, ru1), id2, lp) bud, rud = dw((bu1, ru1), id2, lp)
bdu, rdu = up((bd1, rd1), id2, lp) bdu, rdu = up((bd1, rd1), id2, lp)
...@@ -365,16 +374,14 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, ztw, lp::SpaceParm{4,M,B ...@@ -365,16 +374,14 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, ztw, lp::SpaceParm{4,M,B
gt2 = U[bud,id2,rud] gt2 = U[bud,id2,rud]
end end
if SFBC && (it == 1)
csw[b,ipl,r] = projalg(U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2])) + projalg((U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r]) if !(SFBC && (it == 1))
else
if TWP if TWP
csw[b,ipl,r] = projalg(ztw,U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2])) + projalg(ztw,(U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r]) + csw[b,ipl,r] = projalg(ztw,U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2]) + (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r] +
projalg(ztw,(U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2])) + projalg(ztw,(U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1]) (U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2]) + (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1])
else else
csw[b,ipl,r] = projalg(U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2])) + projalg((U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r]) + csw[b,ipl,r] = projalg(U[b,id1,r]*gt1/(U[b,id2,r]*U[bu2,id1,ru2]) + (U[bd2,id2,rd2]\(U[bd2,id1,rd2]*gt2))/U[b,id1,r] +
projalg((U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2])) + projalg((U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1]) (U[bdd,id2,rdd]*U[bd1,id1,rd1])\(U[bdd,id1,rdd]*U[bd2,id2,rd2]) + (U[b,id2,r]/(U[bd1,id2,rd1]*U[bdu,id1,rdu]))*U[bd1,id1,rd1])
end end
end end
......
src/Groups/GroupSU3.jl
View file @ 51a3f002
...@@ -90,3 +90,4 @@ function isgroup(a::SU3{T}) where T <: AbstractFloat ...@@ -90,3 +90,4 @@ function isgroup(a::SU3{T}) where T <: AbstractFloat
end end
end end
Base.:+(a::SU3{T},b::SU3{T}) where {T} = SU3{T}(a.u11+b.u11,a.u12+b.u12,a.u13+b.u13,a.u21+b.u21,a.u22+b.u22,a.u23+b.u23)
\ No newline at end of file
src/Solvers/Propagators.jl
View file @ 51a3f002
...@@ -2,13 +2,13 @@ ...@@ -2,13 +2,13 @@
""" """
function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, ymws::YMworkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64)
Returns the fermionic progapator for a source at point `y` with color `c` and spin `s`. If the last three arguments are replaced by `time::Int64`, the source is replaced Returns the fermionic progapator for a source at point `y` with color `c` and spin `s`. If the last three arguments are replaced by `time::Int64`, the source is replaced
by a random source in spin and color at t = `time`. by a random source in spin and color at t = `time`.
""" """
function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T} function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, ymws::YMworkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, y::NTuple{4,Int64}, c::Int64, s::Int64) where {T}
function krnlg5!(src) function krnlg5!(src)
b=Int64(CUDA.threadIdx().x) b=Int64(CUDA.threadIdx().x)
...@@ -25,14 +25,14 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, m ...@@ -25,14 +25,14 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, m
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
end end
g5Dw!(pro,U,src,dpar,lp) g5Dw!(pro,U,src,dpar,lp,ymws)
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) CG!(pro,U,DwdagDw!,dpar,lp,dws,ymws,maxiter,tol)
return pro return pro
end end
function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T} function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,ymws::YMworkspace, lp::SpaceParm, maxiter::Int64, tol::Float64, time::Int64) where {T}
function krnlg5!(src) function krnlg5!(src)
b=Int64(CUDA.threadIdx().x) b=Int64(CUDA.threadIdx().x)
...@@ -49,25 +49,25 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, m ...@@ -49,25 +49,25 @@ function propagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm, m
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
end end
g5Dw!(pro,U,src,dpar,lp) g5Dw!(pro,U,src,dpar,lp,ymws)
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
end end
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) CG!(pro,U,DwdagDw!,dpar,lp,dws,ymws,maxiter,tol)
return pro return pro
end end
""" """
function bndpropagator(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, ymws::YMworkspace, 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'.
""" """
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} function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, ymws::YMworkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
function krnlg5!(src) function krnlg5!(src)
b=Int64(CUDA.threadIdx().x) b=Int64(CUDA.threadIdx().x)
...@@ -99,9 +99,9 @@ function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm ...@@ -99,9 +99,9 @@ function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
end end
g5Dw!(pro,U,src,dpar,lp) g5Dw!(pro,U,src,dpar,lp,ymws)
CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) CG!(pro,U,DwdagDw!,dpar,lp,dws,ymws,maxiter,tol)
return pro return pro
end end
test/dirac/test_fp_fa.jl
View file @ 51a3f002
...@@ -4,7 +4,7 @@ using TimerOutputs ...@@ -4,7 +4,7 @@ using TimerOutputs
@timeit "fA_fP test" begin @timeit "fA_fP test" begin
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} function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm{4,6,1,D},ymws, maxiter::Int64, tol::Float64, c::Int64, s::Int64) where {T,D}
function krnlg5!(src) function krnlg5!(src)
b=Int64(CUDA.threadIdx().x) b=Int64(CUDA.threadIdx().x)
...@@ -36,9 +36,9 @@ function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm ...@@ -36,9 +36,9 @@ function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace,lp::SpaceParm
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(src)
end end
g5Dw!(pro,U,src,dpar,lp) g5Dw!(pro,U,src,dpar,lp,ymws)
CG!(pro,U,DwdagDw!,dpar,lp,dws,"ymws","gp",maxiter,tol) CG!(pro,U,DwdagDw!,dpar,lp,dws,ymws,maxiter,tol)
return pro return pro
end end
...@@ -52,6 +52,7 @@ exptheta = exp.(im.*theta./lp.iL); ...@@ -52,6 +52,7 @@ exptheta = exp.(im.*theta./lp.iL);
dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0); dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0);
dws = DiracWorkspace(SU3fund{Float64},Float64,lp); dws = DiracWorkspace(SU3fund{Float64},Float64,lp);
ymws = YMworkspace(SU3, Float64, lp)
U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}));
psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
...@@ -106,6 +107,7 @@ function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1 ...@@ -106,6 +107,7 @@ function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1
dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0); dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0);
dws = DiracWorkspace(SU3fund{Float64},Float64,lp); dws = DiracWorkspace(SU3fund{Float64},Float64,lp);
ymws = YMworkspace(SU3, Float64, lp)
U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}));
psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); psi = scalar_field(Spinor{4,SU3fund{Float64}},lp);
...@@ -113,7 +115,7 @@ function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1 ...@@ -113,7 +115,7 @@ function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1
res = im*zeros(lp.iL[4]) res = im*zeros(lp.iL[4])
for s in 1:4 for c in 1:3 for s in 1:4 for c in 1:3
psi = bndpropagator(U,dpar,dws,lp,1000,prec,c,s); psi = bndpropagator(U,dpar,dws,lp,ymws,1000,prec,c,s);
for t in 1:lp.iL[4] for t in 1:lp.iL[4]
#for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
......
test/dirac/test_solver_plw.jl
View file @ 51a3f002
...@@ -79,8 +79,8 @@ end ...@@ -79,8 +79,8 @@ end
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(pwave) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(pwave)
end end
g5Dw!(prop,U,pwave,dpar,lp) g5Dw!(prop,U,pwave,dpar,lp,ymws)
CG!(prop,U,DwdagDw!,dpar,lp,dws,ymws,gp,10000,1.0e-14) CG!(prop,U,DwdagDw!,dpar,lp,dws,ymws,10000,1.0e-14)
end end
dif = sum(norm2.(prop - prop_th)) dif = sum(norm2.(prop - prop_th))
......
test/dirac/test_solver_rand.jl
View file @ 51a3f002
...@@ -47,11 +47,9 @@ ymws = YMworkspace(SU3, Float64, lp) ...@@ -47,11 +47,9 @@ ymws = YMworkspace(SU3, Float64, lp)
dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0) dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0)
dws = DiracWorkspace(SU3fund{Float64},Float64,lp); dws = DiracWorkspace(SU3fund{Float64},Float64,lp);
randomize!(ymws.mom, lp, ymws) randomize!(ymws.mom, lp, ymws)
U = exp.(ymws.mom) U = exp.(ymws.mom)
rpsi = scalar_field(Spinor{4,SU3fund{Float64}},lp) rpsi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
pfrandomize!(rpsi,lp,dpar) pfrandomize!(rpsi,lp,dpar)
...@@ -69,10 +67,10 @@ end ...@@ -69,10 +67,10 @@ end
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi)
end end
g5Dw!(prop,U,rpsi,dpar,lp) g5Dw!(prop,U,rpsi,dpar,lp,ymws)
CG!(prop,U,DwdagDw!,dpar,lp,dws,ymws,gp,10000,1.0e-14) CG!(prop,U,DwdagDw!,dpar,lp,dws,ymws,10000,1.0e-14)
Dw!(dws.sp,U,prop,dpar,lp) Dw!(dws.sp,U,prop,dpar,lp,ymws)
CUDA.@sync begin CUDA.@sync begin
CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi)
......
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