Fixing the topological charge. Wilson Flow tested, ok

lib/qcd/smearing/WilsonFlow.h

@@ -36,8 +36,10 @@ namespace QCD {
 template <class Gimpl>
 class WilsonFlow: public Smear<Gimpl>{
     unsigned int Nstep;
+    unsigned int measure_interval;
     RealD epsilon;
 
+
     mutable WilsonGaugeAction<Gimpl> SG;
 
     void evolve_step(typename Gimpl::GaugeField&) const;
@@ -47,9 +49,10 @@ class WilsonFlow: public Smear<Gimpl>{
  public:
     INHERIT_GIMPL_TYPES(Gimpl)
 
-    explicit WilsonFlow(unsigned int Nstep, RealD epsilon):
+    explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
         Nstep(Nstep),
         epsilon(epsilon),
+        measure_interval(interval),
         SG(WilsonGaugeAction<Gimpl>(3.0)) {
             // WilsonGaugeAction with beta 3.0
             assert(epsilon > 0.0);
@@ -107,11 +110,20 @@ RealD WilsonFlow<Gimpl>::energyDensityPlaquette(unsigned int step, const GaugeFi
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     out = in;
-    for (unsigned int step = 0; step < Nstep; step++) {
+    for (unsigned int step = 1; step <= Nstep; step++) {
+        auto start = std::chrono::high_resolution_clock::now();
         evolve_step(out);
+        auto end = std::chrono::high_resolution_clock::now();
+        std::chrono::duration<double> diff = end - start;
+        std::cout << "Time to evolve " << diff.count() << " s\n";
         std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
-            << step << " "
+            << step << "  "
             << energyDensityPlaquette(step,out) << std::endl;
+         if( step % measure_interval == 0){
+         std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
+            << step << "  " 
+            << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
+        }
     }
 }
 

lib/qcd/utils/WilsonLoops.h

@@ -197,10 +197,13 @@ public:
 
     std::vector<GaugeMat> U(Nd, grid);
     for (int d = 0; d < Nd; d++) {
+      // this operation is taking too much time
       U[d] = PeekIndex<LorentzIndex>(Umu, d);
     }
     staple = zero;
-    GaugeMat tmp(grid);
+    GaugeMat tmp1(grid);
+    GaugeMat tmp2(grid);
+
 
     for (int nu = 0; nu < Nd; nu++) {
 
@@ -214,22 +217,34 @@ public:
         //      |
         //    __|
         //
+        tmp1 = Cshift(U[nu], mu, 1);
+        tmp2 = Cshift(U[mu], nu, 1);
 
+        staple += tmp1*adj(U[nu] * tmp2);
+
+
+/*
         staple += Gimpl::ShiftStaple(
             Gimpl::CovShiftForward(
                 U[nu], nu,
                 Gimpl::CovShiftBackward(
                     U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
             mu);
-
+*/
         //  __
         // |
         // |__
         //
         //
+
+        tmp2 = adj(U[mu]*tmp1)*U[nu];
+        staple += Cshift(tmp2, nu, -1);
+
+ /*       
         staple += Gimpl::ShiftStaple(
             Gimpl::CovShiftBackward(U[nu], nu,
                                     Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+*/
       }
     }
   }
@@ -289,8 +304,7 @@ public:
       //
       staple = Gimpl::ShiftStaple(
           Gimpl::CovShiftBackward(U[nu], nu,
-                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
+                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-          mu);
     }
   }
 
@@ -307,10 +321,10 @@ public:
       GaugeMat Vup(Umu._grid), Vdn(Umu._grid);
       StapleUpper(Vup, Umu, mu, nu);
       StapleLower(Vdn, Umu, mu, nu);
-      GaugeMat v = adj(Vup) - adj(Vdn);
+      GaugeMat v = Vup - Vdn;
       GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu);  // some redundant copies
       GaugeMat vu = v*u;
-      FS = 0.25*Ta(u*v + Cshift(vu, mu, +1));
+      FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
   }
 
   static Real TopologicalCharge(GaugeLorentz &U){

tests/smearing/Test_WilsonFlow.cc

@@ -42,22 +42,29 @@ int main(int argc, char **argv) {
   GridRedBlackCartesian     RBGrid(latt_size, simd_layout, mpi_layout);
 
   std::vector<int> seeds({1, 2, 3, 4, 5});
+  GridSerialRNG sRNG;
   GridParallelRNG pRNG(&Grid);
   pRNG.SeedFixedIntegers(seeds);
 
   LatticeGaugeField Umu(&Grid), Uflow(&Grid);
   SU<Nc>::HotConfiguration(pRNG, Umu);
+  CheckpointerParameters CPPar("ckpoint_lat", "ckpoint_rng");
+  BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
+
+  CPBin.CheckpointRestore(3000, Umu, sRNG, pRNG);
 
   std::cout << std::setprecision(15);
   std::cout << GridLogMessage << "Plaquette: "
     << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
 
-  WilsonFlow<PeriodicGimplR> WF(200, 0.01);
+  WilsonFlow<PeriodicGimplR> WF(200, 0.01, 50);
 
   WF.smear(Uflow, Umu);
 
   RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
-  std::cout << GridLogMessage << "Plaquette: "<< WFlow_plaq << std::endl;
+  RealD WFlow_TC   = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
+  std::cout << GridLogMessage << "Plaquette         : "<< WFlow_plaq << std::endl;
+  std::cout << GridLogMessage << "TopologicalCharge : "<< WFlow_TC   << std::endl;
 
   std::cout<< GridLogMessage << " Admissibility check:\n";
   const double sp_adm = 0.067;                // admissible threshold