Merge branch 'master' of igit.ific.uv.es:fernando.p.csic.es/latticegpu.jl into fix/flow_obc

src/YM/YM.jl

@@ -165,7 +165,7 @@ include("YMfields.jl")
 export randomize!, zero!, norm2
 
 include("YMact.jl")
-export krnl_plaq!, force_gauge, force_wilson
+export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson
 
 include("YMhmc.jl")
 export gauge_action, hamiltonian, plaquette, HMC!, MD!

src/YM/YMact.jl

@@ -320,6 +320,22 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG,
     return nothing
 end
 
+function bnd_rescale_flw!(frc1, lp::SpaceParm{N,M,BC_OPEN,D}) where {N,M,D}
+
+    @inbounds begin
+        b = Int64(CUDA.threadIdx().x)
+        r = Int64(CUDA.blockIdx().x)
+        I = point_coord((b,r), lp)
+        it = I[N]
+
+        for id in 1:N-1
+            if (((it == 1) || (it == lp.iL[4])))
+                frc1[b,id,r] = 2*frc1[b,id,r]
+            end
+        end
+    end
+    return nothing
+end
 
 ##
 ## SF
@@ -874,7 +890,6 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG,
 end
 
 
-
 ##
 ## PERIODIC
 ##
@@ -1143,6 +1158,38 @@ end
 force_gauge(ymws::YMworkspace, U, c0, gp, lp) = force_gauge(ymws, U, c0, gp.cG[1], gp, lp)
 force_gauge(ymws::YMworkspace, U, gp, lp) = force_gauge(ymws, U, gp.c0, gp.cG[1], gp, lp)
 
+"""
+    function force_gauge_flw(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D})
+
+Computes the force for the gauge flow with Open Boundaries. An aditional factor two in the boundaries
+is included, see
+
+M. Luescher, S. Schaefer: "Lattice QCD with open boundary conditions and twisted-mass reweighting", Comput.Phys.Commun. 184 (2013) 519,
+
+for more details.
+
+"""
+function force_gauge_flw(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}) where {NI,N,M,D}
+
+    ztw = ztwist(gp, lp)
+    if abs(c0-1) < 1.0E-10
+        @timeit "Wilson gauge force" begin
+            force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm)
+        end
+    else
+        @timeit "Improved gauge force" begin
+            force_pln!(ymws.frc1, ymws.frc2, U, gp.Ubnd, cG, ztw, lp::SpaceParm, c0)
+        end
+    end
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz bnd_rescale_flw!(ymws.frc1,lp::SpaceParm)
+    end
+
+    return nothing
+end
+
+
 """ 
     function force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm)
 

src/YM/YMflow.jl

@@ -93,7 +93,7 @@ function Base.show(io::IO, int::FlowIntr{N,T}) where {N,T}
     if N == 0
         println(io, " * Euler schem3")
     elseif N == 1
-        println(io, " * One stage scheme. Coefficients3")
+        println(io, " * One stage scheme. Coefficients")
         println(io, "    stg 1: ", int.e0[1], " ", int.e1[1])
     elseif N == 2
         println(io, " * Two stage scheme. Coefficients:")
@@ -201,6 +201,31 @@ function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpacePar
 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)
 
+function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D}
+    @timeit "Integrating flow equations" begin
+        for i in 1:ns
+            force_gauge_flw(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, 2*eps*int.r)
+
+            for k in 1:NI
+                force_gauge_flw(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
+        end
+    end
+
+    return nothing
+end
+flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D} = flw(U, int, ns, int.eps, gp, lp, ymws)
+
 
 ##
 # Adaptive step size integrators
@@ -320,30 +345,30 @@ Eoft_plaq(U, gp::GaugeParm{T,G,NN}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) w
 
 
 function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd, ztw, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
-    
+
     @inbounds begin
         b = Int64(CUDA.threadIdx().x)
         r = Int64(CUDA.blockIdx().x)
         I = point_coord((b,r), lp)
-        
+
         id1, id2 = lp.plidx[ipl]
         SFBC = ((B == BC_SF_AFWB) || (B == BC_SF_ORBI)) && (id1 == N)
         TWP  = ((I[id1]==1)&&(I[id2]==1))
-        
+
         bu1, ru1 = up((b, r), id1, lp)
         bu2, ru2 = up((b, r), id2, lp)
-        
-        if SFBC && (ru1 != r)
+
+        if SFBC && (point_time((b,r),lp) == lp.iL[end])
             gt = Ubnd[id2]
         else
             gt = U[bu1,id2,ru1]
         end
-        
+
         if TWP
             plx[I] = ztw*tr(U[b,id1,r]*gt / (U[b,id2,r]*U[bu2,id1,ru2]))
         else
             plx[I] = tr(U[b,id1,r]*gt / (U[b,id2,r]*U[bu2,id1,ru2]))
-        end            
+        end
     end
     return nothing
 end

src/YM/YMsf.jl

@@ -92,7 +92,7 @@ end
 """
         function setbndfield(U, phi, lp::SpaceParm)
 
-Sets abelian boundary fields with phases `phi[1]` and `phi[2]` to the configuration `U` at time salice ``x_0=0``.
+Sets abelian boundary fields with phases `phi[1]` and `phi[2]` to the configuration `U` at time slice ``x_0=0``.
 """
 function setbndfield(U, phi, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}