Fermion flow, take one

src/Dirac/Dirac.jl

@@ -480,4 +480,7 @@ export Dw!, DwdagDw!, g5Dw!, pfrandomize!, Csw!
 include("DiracIO.jl")
 export read_prop, save_prop, read_dpar
 
+include("Diracflw.jl")
+export flw
+
 end

src/Dirac/Diracflw.jl

+
+
+import ..YM.flw
+
+function flw(U, psi, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace. dws::DiracWorkspace) where {NI,T}    
+    @timeit "Integrating flow equations" begin
+        for i in 1:ns
+            force_gauge(ymws, U, int.c0, 1, gp, lp)
+
+            DwdagDw!(dws.sAp, U, psi, dpar, dws, dws.st, lp)
+            if int.add_zth
+                add_zth_term(ymws::YMworkspace, U, lp)
+            end
+
+            psi .= psi + 0.5*int.r*eps*dws.sAp
+            ymws.mom .= ymws.frc1
+            U .= expm.(U, ymws.mom, 2*eps*int.r)
+
+            for k in 1:NI
+
+                force_gauge(ymws, U, int.c0, 1, gp, lp)
+                DwdagDw!(dws.sp, U, psi, dpar, dws, dws.st, lp)
+
+                if int.add_zth
+                    add_zth_term(ymws::YMworkspace, U, lp)
+                end
+
+                psi .= psi + 2*eps*int.e0[k].*dws.sAp .+ 2*eps*int.e1[k].*dws.sp
+                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, psi, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace. dws::DiracWorkspace) where {NI,T} = flw(U, psi, int::FlowIntr{NI,T}, ns::Int64, int.eps, gp::GaugeParm, dpar::DiracParam, lp::SpaceParm, ymws::YMworkspace. dws::DiracWorkspace)
+

src/LatticeGPU.jl

@@ -59,6 +59,7 @@ using .Dirac
 export DiracWorkspace, DiracParam
 export Dw!, DwdagDw!, pfrandomize!, g5Dw!, Csw!
 export read_prop, save_prop, read_dpar
+export flw
 
 include("Solvers/Solvers.jl")
 using .Solvers