Added test for the solver

test/dirac/test_solver_plw.jl

+using LatticeGPU, CUDA
+
+#Test for the relation Dw(n|m)^{-1} e^(ipm) = D(p)^{-1} e^{ipn} with a given momenta (if p=0 its randomized), spin and color
+
+function Dwpw_test(;p=0,s=1,c=1)
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
+gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
+ymws = YMworkspace(SU3, Float64, lp)
+dpar = DiracParam{Float64}(SU3fund,1.3,1.0,(1.0,1.0,1.0,1.0),0.0)
+dws = DiracWorkspace(SU3fund{Float64},Float64,lp);
+
+p==0 ? p = Int.(round.(lp.iL.*rand(4),RoundUp)) : nothing
+U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64}))
+
+rm = 2* pi* p./(lp.iL)
+rmom=(rm[1],rm[2],rm[3],rm[4])
+
+pwave = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
+prop = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+prop_th = fill!(scalar_field(Spinor{4,SU3fund{Float64}},lp),zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp))))
+
+#Generate plane wave
+
+for x in 1:lp.iL[1] for y in 1:lp.iL[2] for z in 1:lp.iL[3] for t in 1:lp.iL[4]
+CUDA.@allowscalar pwave[point_index(CartesianIndex{lp.ndim}((x,y,z,t)),lp)...] = exp(im * (x*rmom[1] + y*rmom[2] + z*rmom[3] + t*rmom[4]))*Spinor{4,SU3fund{Float64}}(ntuple(i -> (i==s)*SU3fund{Float64}(ntuple(j -> (j==c)*1.0,3)...),4))
+end end end end
+
+#Th solution
+
+if s == 1
+    vals = (dpar.m0 + 4.0 - sum(cos.(rmom)),0.0,im*sin(rmom[4])+sin(rmom[3]),im*sin(rmom[2])+sin(rmom[1]))
+    for x in 1:lp.iL[1] for y in 1:lp.iL[2] for z in 1:lp.iL[3] for t in 1:lp.iL[4]
+    CUDA.@allowscalar prop_th[point_index(CartesianIndex{lp.ndim}((x,y,z,t)),lp)...] = exp(im * (x*rmom[1] + y*rmom[2] + z*rmom[3] + t*rmom[4]))*
+    ( Spinor{4,SU3fund{Float64}}(ntuple(i -> SU3fund{Float64}(ntuple(j -> (j==c)*vals[i],3)...),4)) )/(sum((sin.(rmom)) .^2) + (dpar.m0+ 4.0 - sum(cos.(rmom)))^2)
+    end end end end
+elseif s == 2
+    vals = (0.0,dpar.m0 + 4.0 - sum(cos.(rmom)),sin(rmom[1]) - im *sin(rmom[2]),-sin(rmom[3])+im*sin(rmom[4]))
+    for x in 1:lp.iL[1] for y in 1:lp.iL[2] for z in 1:lp.iL[3] for t in 1:lp.iL[4]
+    CUDA.@allowscalar prop_th[point_index(CartesianIndex{lp.ndim}((x,y,z,t)),lp)...] = exp(im * (x*rmom[1] + y*rmom[2] + z*rmom[3] + t*rmom[4]))*
+    ( Spinor{4,SU3fund{Float64}}(ntuple(i -> SU3fund{Float64}(ntuple(j -> (j==c)*vals[i],3)...),4)) )/(sum((sin.(rmom)) .^2) + (dpar.m0+ 4.0 - sum(cos.(rmom)))^2)
+    end end end end
+elseif s == 3
+    vals = (-sin(rmom[3])+im*sin(rmom[4]),-sin(rmom[1])-im*sin(rmom[2]),dpar.m0 + 4.0 - sum(cos.(rmom)),0.0)
+    for x in 1:lp.iL[1] for y in 1:lp.iL[2] for z in 1:lp.iL[3] for t in 1:lp.iL[4]
+    CUDA.@allowscalar prop_th[point_index(CartesianIndex{lp.ndim}((x,y,z,t)),lp)...] = exp(im * (x*rmom[1] + y*rmom[2] + z*rmom[3] + t*rmom[4]))*
+    ( Spinor{4,SU3fund{Float64}}(ntuple(i -> SU3fund{Float64}(ntuple(j -> (j==c)*vals[i],3)...),4)) )/(sum((sin.(rmom)) .^2) + (dpar.m0+ 4.0 - sum(cos.(rmom)))^2)
+    end end end end
+else 
+    vals = (-sin(rmom[1])+im*sin(rmom[2]),sin(rmom[3])+im*sin(rmom[4]),0.0,dpar.m0 + 4.0 - sum(cos.(rmom)))
+    for x in 1:lp.iL[1] for y in 1:lp.iL[2] for z in 1:lp.iL[3] for t in 1:lp.iL[4]
+    CUDA.@allowscalar prop_th[point_index(CartesianIndex{lp.ndim}((x,y,z,t)),lp)...] = exp(im * (x*rmom[1] + y*rmom[2] + z*rmom[3] + t*rmom[4]))*
+    ( Spinor{4,SU3fund{Float64}}(ntuple(i -> SU3fund{Float64}(ntuple(j -> (j==c)*vals[i],3)...),4)) )/(sum((sin.(rmom)) .^2) + (dpar.m0+ 4.0 - sum(cos.(rmom)))^2)
+    end end end end
+end
+
+#compute Sum{x} D^-1(x|y)P(y)
+
+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
+
+CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(pwave)
+end
+g5Dw!(prop,U,pwave,dpar,lp)
+CG!(prop,U,DwdagDw!,dpar,lp,dws,10000,1.0e-14)
+
+
+diff = sum(norm2.(prop - prop_th))
+
+if diff < 1.0e-15 
+    print("Test Dplw for s=",s,", c=",c," passed with ",diff,"% error!\n")
+    
+else
+        error("Test Dplw failed with difference: ",diff,"\n")
+end
+
+return nothing
+end
+
+for i in 1:3 for j in 1:4
+    Dwpw_test(c=i,s=j)
+end end
\ No newline at end of file

test/dirac/test_solver_rand.jl

+using CUDA,LatticeGPU
+
+#Check that Dw ( (DwdagDw)^{-1} g5 Dw g5 ) psi = psi for random fields
+
+function pfrandomize!(f,lp::SpaceParm,dpar::DiracParam)
+
+    if dpar.rep == SU3fund && lp.ndim == 4
+        begin
+            p = ntuple(i->CUDA.randn(Float64, lp.bsz, 3, lp.rsz,2),4)
+            CUDA.@sync begin
+                CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_pf!(f,p,lp)
+            end
+        end
+    end
+
+    return nothing
+
+end
+function krnl_assign_pf!(f::AbstractArray{T}, p ,lp::SpaceParm) where {T} 
+    
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+
+            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
+    return nothing
+
+end
+
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0))
+gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0)
+ymws = YMworkspace(SU3, Float64, lp)
+dpar = DiracParam{Float64}(SU3fund,2.3,1.0,(1.0,1.0,1.0,1.0),0.0)
+dws = DiracWorkspace(SU3fund{Float64},Float64,lp);
+
+
+randomize!(ymws.mom, lp, ymws)
+U = exp.(ymws.mom)
+
+
+rpsi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+pfrandomize!(rpsi,lp,dpar)
+
+prop = scalar_field(Spinor{4,SU3fund{Float64}},lp)
+
+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
+
+CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi)
+end
+g5Dw!(prop,U,rpsi,dpar,lp)
+CG!(prop,U,DwdagDw!,dpar,lp,dws,10000,1.0e-14)
+
+Dw!(dws.sp,U,prop,dpar,lp)
+
+CUDA.@sync begin
+    CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi)
+end
+
+res = sum(norm2.(rpsi-dws.sp))
+
+
+if res < 1.0e-15 
+    print("Test Drand passed with ",res,"% error!\n")
+    
+else
+    error("Test Drand failed with difference: ",res,"\n")
+end
+

test/runtests.jl

 
 #include("SAD/test_sad.jl")
 include("flow/test_adapt.jl")
+include("dirac/test_solver_plw.jl")
+include("dirac/test_solver_rand.pl")