From 266ec12b98bb100622ff93ead17c34d2a23f5d0e Mon Sep 17 00:00:00 2001
From: Xiao-Yong Jin <xjin@anl.gov>
Date: Thu, 18 Jul 2024 16:49:07 -0500
Subject: [PATCH] experimental/graph: add stout smear deriv test

---
 src/experimental/graph/gauge.nim   |  7 +++
 src/experimental/graph/tggauge.nim |  2 +-
 src/experimental/graph/tgstout.nim | 98 ++++++++++++++++++++++++++++++
 3 files changed, 106 insertions(+), 1 deletion(-)
 create mode 100644 src/experimental/graph/tgstout.nim

diff --git a/src/experimental/graph/gauge.nim b/src/experimental/graph/gauge.nim
index f058cd6..c1d04ff 100644
--- a/src/experimental/graph/gauge.nim
+++ b/src/experimental/graph/gauge.nim
@@ -15,6 +15,13 @@ type Ggauge* {.final.} = ref object of Gvalue
 
 proc getgauge*(x: Gvalue): Gauge = Ggauge(x).gval
 
+proc update*(x: Gvalue, g: Gauge) =
+  let u = x.getgauge
+  threads:
+    for mu in 0..<u.len:
+      u[mu] := g[mu]
+  x.updated
+
 proc toGvalue*(x: Gauge): Ggauge =
   # proc instead of converter to avoid converting seq
   result = Ggauge(gval: x)
diff --git a/src/experimental/graph/tggauge.nim b/src/experimental/graph/tggauge.nim
index a51ce76..f033597 100644
--- a/src/experimental/graph/tggauge.nim
+++ b/src/experimental/graph/tggauge.nim
@@ -41,7 +41,7 @@ proc ndiff(zt: Gvalue, t: Gscalar): (float, float) =
   (dzdt, e)
 
 template check(ii: tuple[filename:string, line:int, column:int], ast: string, dzdt, e, gdota: float) =
-  if not almostEqual(gdota, dzdt, unitsInLastPlace = 128*1024):
+  if not almostEqual(gdota, dzdt, unitsInLastPlace = 1024*1024):
     checkpoint(ii.filename & ":" & $ii.line & ":" & $ii.column & ": Check failed: " & ast)
     checkpoint("  ndiff: " & $dzdt & " +/- " & $e)
     checkpoint("  grad: " & $gdota)
diff --git a/src/experimental/graph/tgstout.nim b/src/experimental/graph/tgstout.nim
new file mode 100644
index 0000000..e387f98
--- /dev/null
+++ b/src/experimental/graph/tgstout.nim
@@ -0,0 +1,98 @@
+import qex, algorithms/numdiff, gauge/stoutsmear
+import core, scalar, gauge
+
+qexInit()
+
+letParam:
+  lat = @[12,12,12,24]
+  dt = 0.1
+  eps = 0.01
+  nstep = 3
+  beta = 5.4
+  seed:uint = 1234567891
+
+echoParams()
+echo "rank ", myRank, "/", nRanks
+threads: echo "thread ", threadNum, "/", numThreads
+
+let
+  lo = lat.newLayout
+  vol = lo.physVol
+  gc = GaugeActionCoeffs(plaq: beta)
+  g = lo.newgauge
+  u = lo.newgauge
+
+var
+  r = lo.newRNGField(RngMilc6, seed)
+  ss = lo.newStoutSmear(dt)
+
+for i in 0..3:
+  ss.smear(g, g)
+
+g.random r
+g.echoPlaq
+for i in 1..nstep:
+  if i==1:
+    ss.smear(g, u)
+  else:
+    ss.smear(u, u)
+u.echoPlaq
+let sgs = gc.gaugeAction1 u
+echo "smear S: ",sgs
+
+proc act(t: float): float =
+  var ss = lo.newStoutSmear(t)
+  for i in 1..nstep:
+    if i==1:
+      ss.smear(g, u)
+    else:
+      ss.smear(u, u)
+  gc.gaugeAction1 u
+
+var ndt, err: float
+ndiff(ndt, err, act, dt, eps, ordMax=4)
+echo "numdiff smear dS/dt: ",ndt," +/- ",err
+
+proc stout(g, t: Gvalue, n: int): Gvalue =
+  var g = g
+  for i in 1..n:
+    g = axexpmuly(t, gaugeForce(actWilson(-3.0), g), g)
+  g
+
+let
+  gdt = toGvalue dt
+  gg = toGvalue g
+  gs = gg.stout(gdt, nstep)
+  s = gc.gaugeAction gs
+  ddt = s.grad gdt
+
+# echo ddt.treeRepr
+
+gs.eval.getgauge.echoPlaq
+let sgg = s.eval.getfloat
+echo "graph S: ",sgg
+
+let gddt = ddt.eval.getfloat
+echo "graph dS/dt: ",gddt
+# echo ddt.treeRepr
+
+proc gact(t: float): float =
+  gdt.update t
+  s.eval.getfloat
+var gndt, gerr: float
+ndiff(gndt, gerr, gact, dt, eps, ordMax=4)
+echo "numdiff graph dS/dt: ",gndt," +/- ",gerr
+
+let
+  rds = abs((sgs-sgg)/(sgs+sgg))
+  rgdt = abs((ndt-gddt)/(ndt+gddt))
+  rndt = abs((ndt-gndt)/(ndt+gndt))
+echo "rel dS: ",rds
+echo "rel graph dS/dt: ",rgdt
+echo "rel ndiff dS/dt: ",rndt
+
+doassert rds < 1e-10
+doassert rgdt < 1e-10
+doassert rndt < 1e-10
+
+qexFinalize()