From 3c8b39c29119285193a92865cc197156bdd71245 Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Fri, 5 Nov 2021 15:32:28 +0100
Subject: [PATCH] Updated JuliaGPU. Added support for boundary fields

---
 Manifest.toml  |  2 +-
 examples/ym.in |  8 +++++++-
 main/ym.jl     | 26 +++++++++++++++++++-------
 3 files changed, 27 insertions(+), 9 deletions(-)

diff --git a/Manifest.toml b/Manifest.toml
index 37178d9..13fed64 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -175,7 +175,7 @@ version = "0.0.11+0"
 
 [[LatticeGPU]]
 deps = ["CUDA", "Random", "TimerOutputs"]
-git-tree-sha1 = "62d7728203bfe2cf4934648f9cc356cf16f69cf3"
+git-tree-sha1 = "d8e0de54b67807478ead3f1fde4de0a0d94f7277"
 repo-rev = "master"
 repo-url = "https://igit.ific.uv.es/alramos/latticegpu.jl"
 uuid = "958c3683-801a-4582-9cfa-2d6e2ae1763b"
diff --git a/examples/ym.in b/examples/ym.in
index 24bee3f..1f3dc60 100644
--- a/examples/ym.in
+++ b/examples/ym.in
@@ -22,7 +22,13 @@ bcs = "PERIODIC" # Choose PERIODIC, SF
 ntwist = [0,0,0,0,0,0]
 beta = 5.96
 c0   = 1.0 # 1.0 is Wilson action, 1.66... is LW, 3.648 is Iwasaki 
-cG   = 0.0 # Only used with SF bcs
+
+# Parameters for SF boundary conditions.
+# These lines are only needed in case of SF boundary 
+# conditions
+cG   = 0.0 
+phi0 = [-1.047197551196598, 0.0]               # Boundary fields for 
+phiT = [-3.141592653589793, 1.047197551196598] # SF boundary conditions
 
 [HMC]
 integrator = "OMF4" # "LEAPFROPG", "OMF2", "OMF4"
diff --git a/main/ym.jl b/main/ym.jl
index 89ba0b8..a7910e1 100644
--- a/main/ym.jl
+++ b/main/ym.jl
@@ -50,11 +50,20 @@ function read_options(fname)
         GRP = SU3
     end
     T  = Float64
-    gp = GaugeParm{T}(GRP{T},
-                      s["Simulation"]["beta"],
-                      s["Simulation"]["c0"],
-                      (s["Simulation"]["cG"],s["Simulation"]["cG"]))
-    
+    if (BCS == BC_SF_AFWB)
+        gp = GaugeParm{T}(GRP{T},
+                          s["Simulation"]["beta"],
+                          s["Simulation"]["c0"],
+                          (s["Simulation"]["cG"],s["Simulation"]["cG"]),
+                          s["Simulation"]["phiT"], lp.iL[1:end-1])
+        phi0 = s["Simulation"]["phi0"]
+    else
+        gp = GaugeParm{T}(GRP{T},
+                          s["Simulation"]["beta"],
+                          s["Simulation"]["c0"])
+        phi0 = nothing
+    end
+        
     # HMC parameters
     if s["HMC"]["integrator"] == "OMF4"
         int = omf4(T, s["HMC"]["eps"], s["HMC"]["ns"])
@@ -115,16 +124,19 @@ function read_options(fname)
     nthm  = s["HMC"]["nthm"]
     nt    = s["HMC"]["ntraj"]
 
-    return gp, lp, int, flwint, nthm, nt, flw_ns, flw_dtr, fb, flog, GRP, T, BCS
+    return gp, phi0, lp, int, flwint, nthm, nt, flw_ns, flw_dtr, fb, flog, GRP, T, BCS
 end
 
 parsed_args = parse_commandline()
 infile = parsed_args["i"]
 
-gp, lp, int, flwint, nthm, nt, flw_ns, flw_dtr, fb, flog, GRP, T, BCS = read_options(infile)
+gp, phi0, lp, int, flwint, nthm, nt, flw_ns, flw_dtr, fb, flog, GRP, T, BCS = read_options(infile)
 
 U = vector_field(GRP{T}, lp)
 fill!(U, one(GRP{T}))
+if phi0 != nothing
+    setbndfield(U, phi0, lp)
+end
 
 ymws = YMworkspace(GRP, T, lp)
 
-- 
2.24.1