Added run. Updated README

README.md

@@ -28,16 +28,33 @@ julia> Ac = CuArray(1000)
 
 1. To run a simulation use:
 ```bash
-[user@cluster]$ julia --project=/path/to/test_latticegpu.jl /path/to/test_latticegpu.jl/src/ym.jl -i <input file>
+[user@cluster]$ julia --project=/path/to/test_latticegpu.jl /path/to/test_latticegpu.jl/main/ym.jl -i <input file>
 ```
 2. An example input file is available in `examples/ym.in`
 
-### Results: analysis and storage
+### Analysis
 
-Data from testing runs should be collected in the `data` folder. Each
-run should be in its own directory. The input file should be included
-in this git repository. The `log.gz` and `bdio` files should be stored
-in ciclope using `dvc`.
+The `BDIO` file of the simulation can be analized with the
+[ana.jl](./main/ana.jl) program. An example execution is:
+```bash
+[user@cluster]$ julia --project=/path/to/test_latticegpu.jl /path/to/test_latticegpu.jl/main/ana.jl -i <bdio file> -o <dir> [-c <target t^2<E(t)>>]
+```
+The directory `<dir>` should exist before running the program. 
+
+This program only performs some basic analysis, and should be cross-checked.
+
+### Storing results of the simulations
+
+Data from testing runs should be collected in the `data` folder. The
+input file should be included in this git repository. The `log.gz` and
+`bdio` files should be stored in ciclope (path:
+`/net/ciclope-backend/mnt/IFTLattice/ALPHA/test_LatticeGPU`) using
+`dvc`. 
+
+To get a list of all simulation results stored in ciclope, run:
+```bash
+[user@cluster]$ dvc list -R --dvc-only https://gitlab.ift.uam-csic.es/alberto/test_latticegpu.jl.git
+```
 
 ## TODO list
 
@@ -49,11 +66,30 @@ everything works as expected).
 
 ### Reproducing numbers in the literature
 
-#### Flow measurements in arXiv1006.4518
+#### Flow measurements in arXiv:1006.4518
 
 - `SU(3)` at $`\beta = 5.96`$ with lattice size $`24^3 \times 48`$,
   periodic boundary conditions:
   - [input file](./data/SU3_PBC_bt5.96_L24x24x24x48/SU3_PBC_bt5.96_L24x24x24x48.in)
   - [analysis](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis)
   - [results](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis/results.txt)
-  
+- Numerical results for t0/a^2 are barely compatible:
+  - $`\langle t0\rangle = 2.7854(62)`$ (arXiv:1006.4518, table 1).
+  - $`\langle t0\rangle = 2.7719(43)`$ 
+	[results](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis/results.txt)
+- The history of the topological charge does not look good.
+
+#### Flow measurements in arXiv:1710.07526 and arXiv:1805.07106
+
+- These are simulations in `SU(2)` at $`\beta = 2.5`$ with lattice
+  size $`24^3 \times 48`$. 
+  - [input file](./data/SU3_PBC_bt5.96_L24x24x24x48/SU3_PBC_bt5.96_L24x24x24x48.in)
+  - [analysis](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis)
+  - [results](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis/results.txt)
+- The plaquette agrees perfectly (assuming that the reference missed a
+  factor two in the definition of the plaquette):
+  - $`\langle P\rangle = 0.651967(2)`$ (arXiv:1710.07526, table VII).
+  - $`\langle P\rangle = 0.6519656(57)`$
+    ([results](./data/SU3_PBC_bt5.96_L24x24x24x48/analysis/results.txt))
+	
+	

data/SU2_PBC_bt2.5_L24x24x24x48/.gitignore

+/SU2_PBC_bt2.5_L24x24x24x48.bdio
+/SU2_PBC_bt2.5_L24x24x24x48.log.gz

data/SU2_PBC_bt2.5_L24x24x24x48/SU2_PBC_bt2.5_L24x24x24x48.bdio.dvc

+outs:
+- md5: f01bdfa9f7fa56d1cda9801056660440
+  size: 13401066
+  path: SU2_PBC_bt2.5_L24x24x24x48.bdio

data/SU2_PBC_bt2.5_L24x24x24x48/SU2_PBC_bt2.5_L24x24x24x48.in

+#
+# Input file to reproduce results of
+#
+#  Reversibility Violation in the Hybrid Monte Carlo Algorithm
+#  Carsten Urbach
+#  arXiv:1710.07526
+#
+
+[Run]
+user = "alberto"
+host = "som5"
+name = "SU2_PBC_bt2.5_L24x24x24x48" 
+seed = 197803
+
+[Simulation]
+group = "SU2" # Choose SU2, SU3
+dimensions = 4
+size = [24,24,24,48]
+blocks = [4,4,4,4]
+bcs = "PERIODIC" # Choose PERIODIC, SF
+ntwist = [0,0,0,0,0,0]
+beta = 2.5
+c0   = 1.0 # 1.0 is Wilson action, 1.66... is LW, 3.648 is Iwasaki 
+
+[HMC]
+integrator = "OMF4" # "LEAPFROPG", "OMF2", "OMF4"
+ns = 20
+eps = 0.1
+nthm = 1000
+ntraj = 10000
+
+[Flow]
+integrator = "RK3" # EULER, RK2, RK3
+eps = 0.01
+ns  = 550
+dtr = 10 # Frequency of flow measurements
+

data/SU2_PBC_bt2.5_L24x24x24x48/SU2_PBC_bt2.5_L24x24x24x48.log.gz.dvc

+outs:
+- md5: 2c0c082edbe84bb52344cab875afb7ab
+  size: 506897
+  path: SU2_PBC_bt2.5_L24x24x24x48.log.gz

data/SU2_PBC_bt2.5_L24x24x24x48/analysis/results.txt

+
+## Results of analysis of file: data/SU2_PBC_bt2.5_L24x24x24x48/SU2_PBC_bt2.5_L24x24x24x48.bdio
+  - Number of thermalization steps: 1000
+  - Number of measurements:         10000
+  - flow measurements for times:    0.0 - 5.5
+## 
+ # Observable <exp(-dH)>
+  - value:  1.0000006999329316 +/- 2.3958457166358336e-6
+  - tauint: 0.22580553473922516 +/- 0.011413434274011135
+ # Observable <P>
+  - value:  1.3039311752349458 +/- 1.1460898331204123e-5
+  - tauint: 3.303330870452952 +/- 0.48028356376494286
+ # Observable t0 plaquette
+  - value:  5.170751286480197 +/- 0.01711105142052885
+  - tauint: 2.2866696595691676 +/- 0.7424695630327979
+ # Observable w0 plaquette
+  - value:  2.278513607605811 +/- 0.00391718065497019
+  - tauint: 2.2910824986316234 +/- 0.7578526725436543
+ # Observable t0 clover
+  - value:  4.402667266711869 +/- 0.01734714147447308
+  - tauint: 1.6778025209007883 +/- 0.49386821851432905
+ # Observable w0 clover
+  - value:  2.0908926776402916 +/- 0.009592965351155903
+  - tauint: 0.49601062695139997 +/- 0.07843868009171369

main/ana.jl

@@ -29,6 +29,14 @@ function parse_commandline()
         arg_type = String
     end
 
+    @add_arg_table s begin
+        "-c"
+        help = "Value of t^2<E(t)> to determine t0,w0 scales"
+        required = false
+        arg_type = Float64
+        default = 0.3
+    end
+    
     return parse_args(s)
 end
 
@@ -65,11 +73,11 @@ function read_record(dt::Vector, fb, uinfo)
     while BDIO_seek!(fb)
         if BDIO_get_uinfo(fb) == uinfo
             BDIO_read(fb, dt)
-            break
+            return true
         end
     end
 
-    return nothing
+    return false
 end
 
 function read_record(dt::Array{Float64, 2}, fb, uinfo)
@@ -84,11 +92,11 @@ function read_record(dt::Array{Float64, 2}, fb, uinfo)
                     dt[j,i] = foo[j]
                 end
             end
-            break
+            return true
         end
     end
 
-    return nothing
+    return false
 end
 
 function ana_record(fb, nmsm, uinfo)
@@ -108,7 +116,14 @@ end
 function ana_flow_record(fb, nflw, flw_ns, eps, uinfo, tsqet::Union{Nothing, Float64} = nothing)
 
     dt = Array{Float64, 2}(undef, flw_ns+1, nflw)
-    read_record(dt, fb, uinfo)
+    if !read_record(dt, fb, uinfo)
+        if tsqet == nothing
+            return nothing
+        else
+            return nothing, nothing, nothing, nothing
+        end
+    end
+    
     for i in 1:size(dt, 1)
         t = (i-1)*eps
         for j in 1:size(dt,2)
@@ -163,13 +178,25 @@ function ana_flow_record(fb, nflw, flw_ns, eps, uinfo, tsqet::Union{Nothing, Flo
         w0 = sqrt((tsqet-A)/B)
         
         
-        return obs, t0, w0
+        return obs, t0, w0, id
     end
 
 
     return nothing
 end
 
+function plot_et(plq, clv, eps, tsqet)
+
+    t = [(i-1)*eps for i in 1:length(plq)]
+    plot(xlabel="t", ylabel="t^2<E(t)>",legend=:topleft)
+    plot!(t, value.(plq), yerr=err.(plq), label="Plaquette", seriestype=:scatter)
+    plot!(t, value.(clv), yerr=err.(clv), label="Clover", seriestype=:scatter)
+    pl = hline!([tsqet],label="",width=2,color=:black)
+
+    return pl
+end
+
+
 function make_plot(obs)
 
     ens = ensembles(obs)
@@ -193,6 +220,7 @@ end
 parsed_args = parse_commandline()
 fname   = parsed_args["i"]
 dout    = parsed_args["o"]
+tsqet   = parsed_args["c"]
 
 fb = BDIO_open(fname, "r")
 nthm, nmsm, nflw, eps, flw_ns = get_msm_info(fb)
@@ -218,35 +246,49 @@ println(fout, "  - tauint: ", taui(obs, "Ensemble"), " +/- ", dtaui(obs, "Ensemb
 pl = make_plot(obs)
 savefig(pl, dout*"/plaq.pdf")
 
-println(fout, " # Observable t0 plaquette")
-obs, t0, w0 = ana_flow_record(fb, nflw, flw_ns, eps, 12, 0.3)
-uwerr(t0)
-println(fout, "  - value:  ", t0)
-println(fout, "  - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
-pl = make_plot(t0)
-savefig(pl, dout*"/t0_pl.pdf")
-
-println(fout, " # Observable w0 plaquette")
-uwerr(w0)
-println(fout, "  - value:  ", w0)
-println(fout, "  - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
-pl = make_plot(w0)
-savefig(pl, dout*"/w0_pl.pdf")
-
-println(fout, " # Observable t0 clover")
-obs, t0, w0 = ana_flow_record(fb, nflw, flw_ns, eps, 13, 0.3)
-uwerr(t0)
-println(fout, "  - value:  ", t0)
-println(fout, "  - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
-pl = make_plot(t0)
-savefig(pl, dout*"/t0_cl.pdf")
-
-println(fout, " # Observable w0 clover")
-uwerr(w0)
-println(fout, "  - value:  ", w0)
-println(fout, "  - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
-pl = make_plot(w0)
-savefig(pl, dout*"/w0_cl.pdf")
+plq_t0, t0, w0, id = ana_flow_record(fb, nflw, flw_ns, eps, 12, tsqet)
+if t0 == nothing
+    println(fout, " # No flow observables on file")
+else
+    println(fout, " # Observable t0 plaquette")
+    uwerr(t0)
+    println(fout, "  - value:  ", t0)
+    println(fout, "  - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
+    pl = make_plot(t0)
+    savefig(pl, dout*"/t0_pl.pdf")
+
+    println(fout, " # Observable w0 plaquette")
+    uwerr(w0)
+    println(fout, "  - value:  ", w0)
+    println(fout, "  - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
+    pl = make_plot(w0)
+    savefig(pl, dout*"/w0_pl.pdf")
+end
+    
 
+clv_t0, t0, w0, id = ana_flow_record(fb, nflw, flw_ns, eps, 13, tsqet)
+if t0 != nothing
+    println(fout, " # Observable t0 clover")
+    uwerr(t0)
+    println(fout, "  - value:  ", t0)
+    println(fout, "  - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
+    pl = make_plot(t0)
+    savefig(pl, dout*"/t0_cl.pdf")
+    
+    println(fout, " # Observable w0 clover")
+    uwerr(w0)
+    println(fout, "  - value:  ", w0)
+    println(fout, "  - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
+    pl = make_plot(w0)
+    savefig(pl, dout*"/w0_cl.pdf")
+    
+    pl = plot_et(plq_t0, clv_t0, eps, tsqet)
+    savefig(pl, dout*"/tsq_et.pdf")
+    
+    obs = ana_flow_record(fb, nflw, flw_ns, eps, 14, nothing)
+    qtop = obs[id]
+    pl = make_plot(qtop)
+    savefig(pl, dout*"/qtop.pdf")
+end
 
 BDIO_close!(fb)