diff --git a/src/Groups/GroupSU2.jl b/src/Groups/GroupSU2.jl index 7174153..3d98b0a 100644 --- a/src/Groups/GroupSU2.jl +++ b/src/Groups/GroupSU2.jl @@ -21,6 +21,7 @@ dag(a::SU2{T}) where T <: AbstractFloat = inverse(a) norm(a::SU2{T}) where T <: AbstractFloat = sqrt(abs2(a.t1) + abs2(a.t2)) norm2(a::SU2{T}) where T <: AbstractFloat = abs2(a.t1) + abs2(a.t2) tr(g::SU2{T}) where T <: AbstractFloat = complex(2.0*real(g.t1), 0.0) +dev_one(g::SU2{T}) where T <: AbstractFloat = sqrt(( abs2(g.t1 - one(T)) + abs2(g.t2))/2) """ function normalize(a::T) where {T <: Group} diff --git a/src/Groups/GroupSU3.jl b/src/Groups/GroupSU3.jl index c468084..9ac8830 100644 --- a/src/Groups/GroupSU3.jl +++ b/src/Groups/GroupSU3.jl @@ -12,6 +12,9 @@ inverse(a::SU3{T}) where T <: AbstractFloat = SU3{T}(conj(a.u11),conj(a.u21),(a.u12*a.u23 - a.u13*a.u22), conj(a.u12),conj(a.u22),(a.u13*a.u21 - a.u11*a.u23)) dag(a::SU3{T}) where T <: AbstractFloat = inverse(a) tr(a::SU3{T}) where T <: AbstractFloat = a.u11+a.u22+conj(a.u11*a.u22 - a.u12*a.u21) +dev_one(g::SU3{T}) where T <: AbstractFloat = sqrt(( abs2(g.u11 - one(T)) + abs2(g.u12) + abs2(g.u13) + abs2(g.u21) + abs2(g.u22 - one(T)) + abs2(g.u23) )/6) + + function Base.:*(a::SU3{T},b::SU3{T}) where T <: AbstractFloat diff --git a/src/Groups/Groups.jl b/src/Groups/Groups.jl index b049cd9..044d0e6 100644 --- a/src/Groups/Groups.jl +++ b/src/Groups/Groups.jl @@ -47,7 +47,7 @@ include("GroupU1.jl") export U1, U1alg -export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat +export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one end # module diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index a86f9f9..f8184ab 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -17,7 +17,7 @@ include("Groups/Groups.jl") using .Groups export Group, Algebra export SU2, SU2alg, SU3, SU3alg, M3x3, M2x2, U1, U1alg -export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat +export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one include("Space/Space.jl") @@ -41,7 +41,9 @@ using .YM export ztwist export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2 export gauge_action, hamiltonian, plaquette, HMC!, OMF4! -export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3, Eoft_clover, Eoft_plaq, Qtop +export Eoft_clover, Eoft_plaq, Qtop +export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 +export flw, flw_adapt export sfcoupling end # module diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 4f11846..8288c0b 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -137,7 +137,9 @@ include("YMhmc.jl") export gauge_action, hamiltonian, plaquette, HMC!, OMF4! include("YMflow.jl") -export wfl_euler, wfl_rk3, zfl_euler, zfl_rk3, Eoft_clover, Eoft_plaq, Qtop +export FlowIntr, flw, flw_adapt +export Eoft_clover, Eoft_plaq, Qtop +export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 include("YMsf.jl") export sfcoupling diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index c777a3f..3235677 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -9,6 +9,68 @@ ### created: Sat Sep 25 08:37:14 2021 ### + +struct FlowIntr{N,T} + r::T + e0::NTuple{N,T} + e1::NTuple{N,T} + + add_zth::Bool + c0::T + + eps::T + tol::T + eps_ini::T + max_eps::T + sft_fac::T +end + +# pre-defined integrators +wfl_euler(::Type{T}, eps::T, tol::T) where T = FlowIntr{0,T}(one(T),(),(),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) +zfl_euler(::Type{T}, eps::T, tol::T) where T = FlowIntr{0,T}(one(T),(),(),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) +wfl_rk2(::Type{T}, eps::T, tol::T) where T = FlowIntr{1,T}(one(T)/2,(-one(T)/2,),(one(T),),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) +zfl_rk2(::Type{T}, eps::T, tol::T) where T = FlowIntr{1,T}(one(T)/2,(-one(T)/2,),(one(T),),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) +wfl_rk3(::Type{T}, eps::T, tol::T) where T = FlowIntr{2,T}(one(T)/4,(-17/36,-one(T)),(8/9,3/4),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) +zfl_rk3(::Type{T}, eps::T, tol::T) where T = FlowIntr{2,T}(one(T)/4,(-17/36,-one(T)),(8/9,3/4),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) + +function Base.show(io::IO, int::FlowIntr{N,T}) where {N,T} + + if (abs(int.c0-1) < 1.0E-10) + println(io, "WILSON flow integrator") + elseif (abs(int.c0-5/3) < 1.0E-10) && int.add_zth + println(io, "ZEUTHEN flow integrator") + elseif (abs(int.c0-5/3) < 1.0E-10) && !int.add_zth + println(io, "SYMANZIK flow integrator") + else + println(io, "CUSTOM flow integrator") + if int.add_zth + println(io, " - ", int.c0, " (with zeuthen term)") + else + println(io, " - ", int.c0) + end + end + + if N == 0 + println(io, " * Euler schem3") + elseif N == 1 + println(io, " * One stage scheme. Coefficients3") + println(io, " stg 1: ", int.e0[1], " ", int.e1[1]) + elseif N == 2 + println(io, " * Two stage scheme. Coefficients:") + println(io, " stg 1: ", int.e0[1], " ", int.e1[1]) + println(io, " stg 2: ", int.e0[2], " ", int.e1[2]) + end + + println(io, " * Fixed step size parameters: eps = ", int.eps) + println(io, " * Adaptive step size parameters: tol = ", int.tol) + println(io, " - max eps: ", int.max_eps) + println(io, " - initial eps: ", int.eps_ini) + println(io, " - safety scale: ", int.sft_fac) + + return nothing +end + + """ function add_zth_term(ymws::YMworkspace, U, lp) @@ -73,69 +135,82 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S return nothing end - -function flw_euler(U, ns, eps, c0, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace; add_zth=false) - - @timeit "Integrating flow equations (Euler)" begin +function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} + @timeit "Integrating flow equations" begin for i in 1:ns - force_gauge(ymws, U, c0, 1, gp, lp) - if add_zth - add_zth_term(ymws::YMworkspace, U, lp) - end - U .= expm.(U, ymws.frc1, 2*eps) - end - end - - return nothing -end - -function flw_rk3(U, ns, eps, c0, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace; add_zth=false) - - @timeit "Integrating flow equations (RK3)" begin - for i in 1:ns - e0 = eps/2 - force_gauge(ymws, U, c0, 1, gp, lp) - if add_zth + force_gauge(ymws, U, int.c0, 1, gp, lp) + if int.add_zth add_zth_term(ymws::YMworkspace, U, lp) end ymws.mom .= ymws.frc1 - U .= expm.(U, ymws.mom, e0) - - e0 = -34*eps/36 - e1 = 16*eps/9 - force_gauge(ymws, U, c0, 1, gp, lp) - if add_zth - add_zth_term(ymws::YMworkspace, U, lp) + U .= expm.(U, ymws.mom, 2*eps*int.r) + + for k in 1:NI + force_gauge(ymws, U, int.c0, 1, gp, lp) + if int.add_zth + add_zth_term(ymws::YMworkspace, U, lp) + end + ymws.mom .= int.e0[k].*ymws.mom .+ int.e1[k].*ymws.frc1 + U .= expm.(U, ymws.mom, 2*eps) end - ymws.mom .= e0.*ymws.mom .+ e1.*ymws.frc1 - U .= expm.(U, ymws.mom) - - e1 = 6*eps/4 - force_gauge(ymws, U, c0, 1, gp, lp) - if add_zth - add_zth_term(ymws::YMworkspace, U, lp) - end - ymws.mom .= e1.*ymws.frc1 .- ymws.mom - U .= expm.(U, ymws.mom) end end - + return nothing end +flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} = flw(U, int, ns, int.eps, gp, lp, ymws) +## +# Adaptive step size integrators +## - -wfl_euler(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 1, gp, lp, ymws) -zfl_euler(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_euler(U, ns, eps, 5.0/3.0, gp, lp, ymws, add_zth=true) -wfl_rk3(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 1, gp, lp, ymws) -zfl_rk3(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 5.0/3.0, gp, lp, ymws, add_zth=true) +function flw_adapt(U, int::FlowIntr{NI,T}, tend::T, epsini::T, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} + + eps = int.eps_ini + dt = tend + nstp = 0 + while true + ns = convert(Int64, floor(dt/eps)) + if ns > 10 + flw(U, int, 9, eps, gp, lp, ymws) + ymws.U1 .= U + flw(U, int, 2, eps/2, gp, lp, ymws) + flw(ymws.U1, int, 1, eps, gp, lp, ymws) + + dt = dt - 10*eps + nstp = nstp + 10 + + # adjust step size + ymws.U1 .= ymws.U1 ./ U + maxd = CUDA.mapreduce(dev_one, max, ymws.U1, init=zero(tend)) + eps = min(int.max_eps, 2*eps, int.sft_fac*eps*(int.tol/maxd)^(one(tend)/3)) + + else + flw(U, int, ns, eps, gp, lp, ymws) + dt = dt - ns*eps + + flw(U, int, 1, dt, gp, lp, ymws) + dt = zero(tend) + + nstp = nstp + ns + 1 + end + + if dt == zero(tend) + break + end + end + + return nstp, eps +end +flw_adapt(U, int::FlowIntr{NI,T}, tend::T, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} = flw_adapt(U, int, tend, int.eps_ini, gp, lp, ymws) ## # Observables ## + """ function Eoft_plaq([Eslc,] U, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) @@ -151,7 +226,7 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMwo tp = ntuple(i->i, N-1) V3 = prod(lp.iL[1:end-1]) - + fill!(Eslc,zero(T)) Etmp = zeros(T,lp.iL[end]) for ipl in 1:M @@ -159,8 +234,8 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMwo CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq_pln!(ymws.cm, U, gp.Ubnd, ztw[ipl], ipl, lp) end - - Etmp .= (gp.ng .- reshape(Array(CUDA.mapreduce(real, +, ymws.cm;dims=tp)),lp.iL[end])/V3 ) + + Etmp .= (gp.ng .- reshape(Array(CUDA.mapreduce(real, +, ymws.cm;dims=tp)),lp.iL[end]) ./ V3 ) if ipl < N for it in 2:lp.iL[end] Eslc[it,ipl] = Etmp[it] + Etmp[it-1] @@ -201,7 +276,7 @@ function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd::T, ztw, ipl, lp::SpacePa end I = point_coord((b,r), lp) - plx[I] = ztw*tr(U[b,1,r]*gt / (U[b,2,r]*U[bu2,id1,ru2])) + plx[I] = ztw*tr(U[b,id1,r]*gt / (U[b,id2,r]*U[bu2,id1,ru2])) return nothing end @@ -219,22 +294,23 @@ function Qtop(Qslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace) ztw = ztwist(gp, lp) tp = (1,2,3) + fill!(ymws.rm, zero(eltype(ymws.rm))) CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,6, ztw[1], ztw[5], lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,6, ztw[1], ztw[6], lp) end CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp) end CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 2,5, ztw[2], ztw[4], lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 2,5, ztw[2], ztw[5], lp) end CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, +, ymws.frc1, ymws.frc2, U, lp) end CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,4, ztw[3], ztw[6], lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,4, ztw[3], ztw[4], lp) end CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp) diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 9602e90..1a3dd8d 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -37,7 +37,7 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) return S end -function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) +function plaquette(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D} ztw = ztwist(gp, lp) @timeit "Plaquette measurement" begin @@ -45,7 +45,7 @@ function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_plaq!(ymws.cm, U, gp.Ubnd, one(gp.cG[1]), ztw, lp) end end - + return CUDA.mapreduce(real, +, ymws.cm)/(prod(lp.iL)*lp.npls) end diff --git a/src/main/times.jl b/src/main/times.jl index ce3069d..7743fab 100644 --- a/src/main/times.jl +++ b/src/main/times.jl @@ -7,7 +7,7 @@ Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU") using LatticeGPU # Set lattice/block size -ntwist = (0,0,0,1,0,0) +ntwist = (0,0,0,0,0,0) lp = SpaceParm{4}((32,32,32,32), (4,4,4,4), BC_PERIODIC, ntwist) println("Space Parameters: ", lp) @@ -44,6 +44,8 @@ println("\n## WILSON ACTION/FLOW TIMES") gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 1.0, (0.5,0.5)) println("Gauge Parameters: ", gp) +flwint = wfl_rk3(PREC, 0.005, 1.0E-6) +println(flwint) println("Initial Action: ") @time S = gauge_action(U, lp, gp, ymws) @@ -62,11 +64,12 @@ for i in 1:4 # println("SF coupling: ", x, " ", y) end -wfl_rk3(U, 1, 0.01, gp, lp, ymws) - +Ucp = Array(U) println("Action: ", gauge_action(U, lp, gp, ymws)) +flw(U, flwint, 1, gp, lp, ymws) + println("Time for 100 steps of RK3 flow integrator: ") -@time wfl_rk3(U, 100, 0.01, gp, lp, ymws) +@time flw(U, flwint, 400, gp, lp, ymws) eoft = Eoft_plaq(U, gp, lp, ymws) eoft = Eoft_clover(U, gp, lp, ymws) qtop = Qtop(U, gp, lp, ymws) @@ -77,8 +80,6 @@ println("Plaq: ", eoft) println("Clov: ", eoft) @time qtop = Qtop(U, gp, lp, ymws) println("Qtop: ", qtop) - -println("Action: ", gauge_action(U, lp, gp, ymws)) println("## END Wilson action/flow measurements") # Set gauge parameters @@ -87,6 +88,9 @@ println("\n## IMPROVED ACTION/FLOW TIMES") gp = GaugeParm{PREC}(GRP{PREC}, 6.0, 5/6, (0.5,0.5)) println("Gauge Parameters: ", gp) +flwint = zfl_rk3(PREC, 0.01, 1.0E-6) +println(flwint) + println("Initial Action: ") @time S = gauge_action(U, lp, gp, ymws) @@ -100,11 +104,11 @@ for i in 1:4 println("# Plaquette: ", pl[end], "\n") end -zfl_rk3(U, 1, 0.01, gp, lp, ymws) +flw(U, flwint, 1, gp, lp, ymws) println("Action: ", gauge_action(U, lp, gp, ymws)) println("Time for 100 steps of RK3 flow integrator: ") -@time zfl_rk3(U, 100, 0.01, gp, lp, ymws) +@time flw(U, flwint, 100, gp, lp, ymws) println("Action: ", gauge_action(U, lp, gp, ymws)) println("## END improved action/flow measurements")