Improved support for flow integrators

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}

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
 

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

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

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

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(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T}    
-function flw_euler(U, ns, eps, c0, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace; add_zth=false)
+    @timeit "Integrating flow equations" begin
-    
-    @timeit "Integrating flow equations (Euler)" begin
         for i in 1:ns
-            force_gauge(ymws, U, c0, 1, gp, lp)
+            force_gauge(ymws, U, int.c0, 1, gp, lp)
-            if add_zth
+            if int.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
                 add_zth_term(ymws::YMworkspace, U, lp)
             end
             ymws.mom .= ymws.frc1
-            U .= expm.(U, ymws.mom, e0)
+            U .= expm.(U, ymws.mom, 2*eps*int.r)
-            
+
-            e0 = -34*eps/36
+            for k in 1:NI
-            e1 = 16*eps/9
+                force_gauge(ymws, U, int.c0, 1, gp, lp)
-            force_gauge(ymws, U, c0, 1, gp, lp)
+                if int.add_zth
-            if add_zth
+                    add_zth_term(ymws::YMworkspace, U, lp)
-                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
+##
 
-                              
+function flw_adapt(U, int::FlowIntr{NI,T}, tend::T, epsini::T, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T}
-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)
+    eps = int.eps_ini
-wfl_rk3(U, ns, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 1, gp, lp, ymws)
+    dt  = tend
-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)
+    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)

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
     

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")