SFORBI different from SFAFWB.

This commit is contained in:
Fernando P.Panadero 2024-07-04 15:22:20 +02:00
parent bc06079664
commit 4c560ab0fe

View file

@ -324,13 +324,99 @@ end
## ##
## SF ## SF
## ##
function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::NTuple{NB,T}, cG, ztw, lp::Union{SpaceParm{N,M,BC_SF_ORBI,D},SpaceParm{N,M,BC_SF_AFWB,D}}) where {T,NB,N,M,D} function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::NTuple{NB,T}, cG, ztw, lp::SpaceParm{N,M,BC_SF_ORBI,D}) where {T,NB,N,M,D}
b = Int64(CUDA.threadIdx().x) b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x) r = Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp) I = point_coord((b,r), lp)
it = I[N] it = I[N]
ipl = 0
S = zero(eltype(plx))
@inbounds begin
for id1 in N:-1:1
bu1, ru1 = up((b, r), id1, lp)
SFBC = (id1==N)
for id2 = 1:id1-1
bu2, ru2 = up((b, r), id2, lp)
ipl = ipl + 1
TWP = (I[id1]==1) && (I[id2]==1)
TWH2 = TWP || ( (I[id1]==lp.iL[id1]) && (I[id2]==1) )
TWH3 = TWP || ( (I[id1]==1) && (I[id2]==lp.iL[id2]) )
# H2 staple
(b1, r1) = up((b,r), id1, lp)
(b2, r2) = up((b1,r1), id1, lp)
if SFBC && (it == lp.iL[end]-1)
gb = Ubnd[id2]
else
gb = U[b2,id2,r2]
end
(b2, r2) = up((b1,r1), id2, lp)
h2 = (U[b1,id1,r1]*gb)/U[b2,id1,r2]
# H3 staple
(b1, r1) = up((b,r), id2, lp)
(b2, r2) = up((b1,r1), id2, lp)
(b3, r3) = up((b1,r1), id1, lp)
if SFBC && (it == lp.iL[end])
gc = Ubnd[id2]
else
gc = U[b3,id2,r3]
end
h3 = (U[b1,id2,r1]*U[b2,id1,r2])/gc
# END staples
if SFBC && (it == lp.iL[end])
ga = Ubnd[id2]
else
ga = U[bu1,id2,ru1]
end
g2 = U[b,id2,r]\U[b,id1,r]
if (it == lp.iL[end]) && SFBC
S += cG*(c0*tr(g2*ga/U[bu2,id1,ru2])) + (c1)*tr(g2*ga/h3) + (c1/2)*tr((g2*ga/U[bu2,id1,ru2])*g2*ga/U[bu2,id1,ru2])
elseif (it == 1) && SFBC
S += cG*(c0*tr(g2*ga/U[bu2,id1,ru2])) + (c1)*tr(g2*ga/h3) + c1*tr(g2*h2/U[bu2,id1,ru2]) + (c1/2)*tr((g2*ga/U[bu2,id1,ru2])*g2*ga/U[bu2,id1,ru2])
else
if TWP
S += (ztw[ipl]*c0)*tr(g2*ga/U[bu2,id1,ru2])
else
S += c0*tr(g2*ga/U[bu2,id1,ru2])
end
if TWH2
S += (ztw[ipl]*c1)*tr(g2*h2/U[bu2,id1,ru2])
else
S += c1*tr(g2*h2/U[bu2,id1,ru2])
end
if TWH3
S += (ztw[ipl]*c1)*tr(g2*ga/h3)
else
S += c1*tr(g2*ga/h3)
end
end
end
end
plx[I] = S
end
return nothing
end
function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, Ubnd::NTuple{NB,T}, cG, ztw, lp::SpaceParm{N,M,BC_SF_AFWB,D}) where {T,NB,N,M,D}
b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
it = I[N]
ipl = 0 ipl = 0
S = zero(eltype(plx)) S = zero(eltype(plx))
@inbounds begin @inbounds begin
@ -509,7 +595,146 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, Ubnd, cG, ztw,
return nothing return nothing
end end
function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, ztw, ipl, lp::Union{SpaceParm{N,M,BC_SF_ORBI,D},SpaceParm{N,M,BC_SF_AFWB,D}}) where {T,N,M,D} function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, ztw, ipl, lp::SpaceParm{N,M,BC_SF_ORBI,D}) where {T,N,M,D}
b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x)
I = point_coord((b,r), lp)
it = I[N]
@inbounds begin
id1, id2 = lp.plidx[ipl]
bu1, ru1 = up((b, r), id1, lp)
bu2, ru2 = up((b, r), id2, lp)
SFBC = (id1 == N)
TWP = (I[id1]==1) && (I[id2]==1)
TWH1 = TWP || ( (I[id1]==1) && (I[id2]==2) )
TWH2 = TWP || ( (I[id1]==lp.iL[id1]) && (I[id2]==1) )
TWH3 = TWP || ( (I[id1]==1) && (I[id2]==lp.iL[id2]) )
TWH4 = TWP || ( (I[id1]==2) && (I[id2]==1) )
# H1 staple
(b1, r1) = dw((b,r), id2, lp)
(b2, r2) = up((b1,r1), id1, lp)
if SFBC && (it == lp.iL[end])
gc = Ubnd[id2]
else
gc = U[b2,id2,r2]
end
h1 = (U[b1,id2,r1]\U[b1,id1,r1])*gc
# H2 staple
(b1, r1) = up((b,r), id1, lp)
(b2, r2) = up((b1,r1), id1, lp)
if SFBC && (it == lp.iL[end]-1)
gb = Ubnd[id2]
else
gb = U[b2,id2,r2]
end
(b2, r2) = up((b1,r1), id2, lp)
h2 = (U[b1,id1,r1]*gb)/U[b2,id1,r2]
# H3 staple
(b1, r1) = up((b,r), id2, lp)
(b2, r2) = up((b1,r1), id2, lp)
(b3, r3) = up((b1,r1), id1, lp)
if SFBC && (it == lp.iL[end])
gc = Ubnd[id2]
else
gc = U[b3,id2,r3]
end
h3 = (U[b1,id2,r1]*U[b2,id1,r2])/gc
# H4 staple
(b1, r1) = dw((b,r), id1, lp)
(b2, r2) = up((b1,r1), id2, lp)
h4 = (U[b1,id1,r1]\U[b1,id2,r1])*U[b2,id1,r2]
# END staples
if SFBC && (it == lp.iL[end])
ga = Ubnd[id2]
else
ga = U[bu1,id2,ru1]
end
g1 = ga/U[bu2,id1,ru2]
g2 = U[b,id2,r]\U[b,id1,r]
if SFBC && (it == 1)
X = (cG*c0)*projalg(U[b,id1,r]*g1/U[b,id2,r]) + c1*projalg(U[b,id1,r]*h2/(U[b,id2,r]*U[bu2,id1,ru2])) +
(c1)*projalg(U[b,id1,r]*ga/(U[b,id2,r]*h3)) + c1*(projalg(U[b,id1,r]*g1*g2*g1/U[b,id2,r]))
frc1[b,id1,r] -= X
frc2[bu1,id2,ru1] -= (cG*c0)*projalg(g1*g2) + c1*projalg((ga/h3)*g2) +
c1*projalg((g1/U[b,id2,r])*h1) + c1*projalg(g1*g2*g1*g2)
frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + c1*projalg((U[b,id2,r]\h1)*g1) +
c1*projalg(g2*h2/U[bu2,id1,ru2]) + c1*projalg(g2*g1*g2*g1)
elseif SFBC && (it == lp.iL[end])
X = cG*c0*projalg(U[b,id1,r]*g1/U[b,id2,r]) +
c1 * (projalg(U[b,id1,r]*ga/(U[b,id2,r]*h3))) + c1*(projalg(U[b,id1,r]*g1*g2*g1/U[b,id2,r]))
frc1[b,id1,r] -= X + c1*projalg(U[b,id1,r]*g1/h4)
frc1[b,id2,r] += X + c1*projalg(h1*g1/U[b,id2,r])
frc2[bu2,id1,ru2] += (cG*c0)*projalg(g2*g1) + c1*projalg((U[b,id2,r]\h1)*g1) +
c1 * projalg(h4\U[b,id1,r]*g1) + c1* projalg(g2*g1*g2*g1)
else
if TWP
X = projalg(c0*ztw,U[b,id1,r]*g1/U[b,id2,r])
frc2[bu1,id2,ru1] -= projalg(c0*ztw,g1*g2)
frc2[bu2,id1,ru2] += projalg(c0*ztw,g2*g1)
else
X = c0*projalg(U[b,id1,r]*g1/U[b,id2,r])
frc2[bu1,id2,ru1] -= c0*projalg(g1*g2)
frc2[bu2,id1,ru2] += c0*projalg(g2*g1)
end
if TWH1
frc1[b,id2,r] += projalg(ztw*c1,h1*g1/U[b,id2,r])
frc2[bu1,id2,ru1] -= projalg(ztw*c1,(g1/U[b,id2,r])*h1)
frc2[bu2,id1,ru2] += projalg(ztw*c1,(U[b,id2,r]\h1)*g1)
else
frc1[b,id2,r] += c1*projalg(h1*g1/U[b,id2,r])
frc2[bu1,id2,ru1] -= c1*projalg((g1/U[b,id2,r])*h1)
frc2[bu2,id1,ru2] += c1*projalg((U[b,id2,r]\h1)*g1)
end
if TWH2
X += projalg(ztw*c1,U[b,id1,r]*h2/(U[b,id2,r]*U[bu2,id1,ru2]))
frc2[bu2,id1,ru2] += projalg(ztw*c1,g2*h2/U[bu2,id1,ru2])
else
X += c1*projalg(U[b,id1,r]*h2/(U[b,id2,r]*U[bu2,id1,ru2]))
frc2[bu2,id1,ru2] += c1*projalg(g2*h2/U[bu2,id1,ru2])
end
if TWH3
X += projalg(ztw*c1,U[b,id1,r]*ga/(U[b,id2,r]*h3))
frc2[bu1,id2,ru1] -= projalg(ztw*c1,(ga/h3)*g2)
else
X += c1*projalg(U[b,id1,r]*ga/(U[b,id2,r]*h3))
frc2[bu1,id2,ru1] -= c1*projalg((ga/h3)*g2)
end
if TWH4
frc1[b,id1,r] -= projalg(ztw*c1,U[b,id1,r]*g1/h4)
frc2[bu1,id2,ru1] -= projalg(ztw*c1,(g1/h4)*U[b,id1,r])
frc2[bu2,id1,ru2] += projalg(ztw*c1,h4\U[b,id1,r]*g1)
else
frc1[b,id1,r] -= c1*projalg(U[b,id1,r]*g1/h4)
frc2[bu1,id2,ru1] -= c1*projalg((g1/h4)*U[b,id1,r])
frc2[bu2,id1,ru2] += c1*projalg(h4\U[b,id1,r]*g1)
end
frc1[b,id1,r] -= X
frc1[b,id2,r] += X
end
end
return nothing
end
function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, ztw, ipl, lp::SpaceParm{N,M,BC_SF_AFWB,D}) where {T,N,M,D}
b = Int64(CUDA.threadIdx().x) b = Int64(CUDA.threadIdx().x)
r = Int64(CUDA.blockIdx().x) r = Int64(CUDA.blockIdx().x)