From aaf1e33a77a31ba0d06e69d1ed476722ce015f5c Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Fri, 2 Jun 2017 16:32:35 +0100
Subject: [PATCH] Adding adaptive integration in the WilsonFlow

---
 lib/qcd/smearing/WilsonFlow.h     | 53 +++++++++++++++++++++++++++++--
 tests/smearing/Test_WilsonFlow.cc | 27 ++++++++++++++++
 2 files changed, 78 insertions(+), 2 deletions(-)

diff --git a/lib/qcd/smearing/WilsonFlow.h b/lib/qcd/smearing/WilsonFlow.h
index 00b7622e..1053775d 100644
--- a/lib/qcd/smearing/WilsonFlow.h
+++ b/lib/qcd/smearing/WilsonFlow.h
@@ -37,15 +37,15 @@ template <class Gimpl>
 class WilsonFlow: public Smear<Gimpl>{
     unsigned int Nstep;
     unsigned int measure_interval;
-    RealD epsilon;
+    mutable RealD epsilon, taus;
 
 
     mutable WilsonGaugeAction<Gimpl> SG;
 
     void evolve_step(typename Gimpl::GaugeField&) const;
+    void evolve_step_adaptive(typename Gimpl::GaugeField&) const;
     RealD tau(unsigned int t)const {return epsilon*(t+1.0); }
 
-
  public:
     INHERIT_GIMPL_TYPES(Gimpl)
 
@@ -101,20 +101,63 @@ void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U) const{
     Gimpl::update_field(Z, U, -2.0*epsilon);    // V(t+e) = exp(ep*Z)*W2
 }
 
+template <class Gimpl>
+void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U) const {
+    GaugeField Z(U._grid);
+    GaugeField Zprime(U._grid);
+    GaugeField tmp(U._grid), Uprime(U._grid);
+    Uprime = U;
+    SG.deriv(U, Z);
+    Zprime = -Z;
+    Z *= 0.25;                                  // Z0 = 1/4 * F(U)
+    Gimpl::update_field(Z, U, -2.0*epsilon);    // U = W1 = exp(ep*Z0)*W0
+
+    Z *= -17.0/8.0;
+    SG.deriv(U, tmp); Z += tmp;                 // -17/32*Z0 +Z1
+    Zprime += 2.0*tmp;
+    Z *= 8.0/9.0;                               // Z = -17/36*Z0 +8/9*Z1
+    Gimpl::update_field(Z, U, -2.0*epsilon);    // U_= W2 = exp(ep*Z)*W1
+    
+
+    Z *= -4.0/3.0;
+    SG.deriv(U, tmp); Z += tmp;                 // 4/3*(17/36*Z0 -8/9*Z1) +Z2
+    Z *= 3.0/4.0;                               // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
+    Gimpl::update_field(Z, U, -2.0*epsilon);    // V(t+e) = exp(ep*Z)*W2
+
+    // Ramos 
+    Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
+    // Compute distance as norm^2 of the difference
+    GaugeField diffU = U - Uprime;
+    RealD diff = norm2(diffU);
+    // adjust integration step
+    
+    taus += epsilon;
+    std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
+    
+    epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
+    std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
+
+}
+
 template <class Gimpl>
 RealD WilsonFlow<Gimpl>::energyDensityPlaquette(unsigned int step, const GaugeField& U) const {
     RealD td = tau(step);
+    //RealD td = tau;
     return 2.0 * td * td * SG.S(U)/U._grid->gSites();
 }
 
 
 //#define WF_TIMING 
 
+
+
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     out = in;
+    taus = epsilon;
     for (unsigned int step = 1; step <= Nstep; step++) {
         auto start = std::chrono::high_resolution_clock::now();
+        std::cout << GridLogMessage << "Evolution time :"<< tau(step) << std::endl;
         evolve_step(out);
         auto end = std::chrono::high_resolution_clock::now();
         std::chrono::duration<double> diff = end - start;
@@ -130,6 +173,12 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
             << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
         }
     }
+    std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
+              << step << "  "
+              << energyDensityPlaquette(Nstep,out) << std::endl;
+    std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
+              << step << "  " 
+              << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
 }
 
 }  // namespace QCD
diff --git a/tests/smearing/Test_WilsonFlow.cc b/tests/smearing/Test_WilsonFlow.cc
index c337bf29..77e7c5d1 100644
--- a/tests/smearing/Test_WilsonFlow.cc
+++ b/tests/smearing/Test_WilsonFlow.cc
@@ -101,7 +101,9 @@ int main(int argc, char **argv) {
 
   RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
   RealD WFlow_TC   = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
+  RealD WFlow_T0   = WF.energyDensityPlaquette(WFPar.steps, Uflow);
   std::cout << GridLogMessage << "Plaquette          "<< conf << "   " << WFlow_plaq << std::endl;
+  std::cout << GridLogMessage << "T0                 "<< conf << "   " << WFlow_T0 << std::endl;
   std::cout << GridLogMessage << "TopologicalCharge  "<< conf << "   " << WFlow_TC   << std::endl;
 
   std::cout<< GridLogMessage << " Admissibility check:\n";
@@ -121,3 +123,28 @@ int main(int argc, char **argv) {
   }
   Grid_finalize();
 }  // main
+
+
+/*
+Input file example
+
+
+JSON
+
+{
+    "WilsonFlow":{
+	"steps": 200,
+	"step_size": 0.01,
+	"meas_interval": 50
+    },
+    "Configurations":{
+	"conf_prefix": "ckpoint_lat",
+	"rng_prefix": "ckpoint_rng",
+	"StartConfiguration": 3000,
+	"EndConfiguration": 3000,
+	"Skip": 5
+    }
+}
+
+
+*/
\ No newline at end of file