Merge branch 'develop' into feature/hadrons

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&, RealD);
     RealD tau(unsigned int t)const {return epsilon*(t+1.0); }
 
-
  public:
     INHERIT_GIMPL_TYPES(Gimpl)
 
@@ -75,7 +75,9 @@ class WilsonFlow: public Smear<Gimpl>{
         // undefined for WilsonFlow
     }
 
+    void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau);
     RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const;
+    RealD energyDensityPlaquette(const GaugeField& U) const;
 };
 
 
@@ -101,20 +103,70 @@ 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, RealD maxTau) {
+    if (maxTau - taus < epsilon){
+        epsilon = maxTau-taus;
+    }
+    std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
+    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);
     return 2.0 * td * td * SG.S(U)/U._grid->gSites();
 }
 
+template <class Gimpl>
+RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const GaugeField& U) const {
+    return 2.0 * taus * taus * SG.S(U)/U._grid->gSites();
+}
+
 
 //#define WF_TIMING 
 
+
+
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     out = in;
     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;
@@ -132,6 +184,30 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     }
 }
 
+template <class Gimpl>
+void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
+    out = in;
+    taus = epsilon;
+    unsigned int step = 0;
+    do{
+        step++;
+        std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
+        evolve_step_adaptive(out, maxTau);
+        std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
+            << step << "  "
+            << energyDensityPlaquette(out) << std::endl;
+         if( step % measure_interval == 0){
+         std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
+            << step << "  " 
+            << WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
+        }
+    } while (taus < maxTau);
+
+
+
+}
+
+
 }  // namespace QCD
 }  // namespace Grid
 

tests/smearing/Test_WilsonFlow.cc

@@ -33,7 +33,8 @@ namespace Grid{
     GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
             int, steps,
             double, step_size,
-            int, meas_interval);
+            int, meas_interval,
+            double, maxTau); // for the adaptive algorithm
        
 
     template <class ReaderClass >
@@ -97,11 +98,13 @@ int main(int argc, char **argv) {
 
   WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
 
-  WF.smear(Uflow, Umu);
+  WF.smear_adaptive(Uflow, Umu, WFPar.maxTau);
 
   RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
   RealD WFlow_TC   = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
+  RealD WFlow_T0   = WF.energyDensityPlaquette(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 +124,29 @@ int main(int argc, char **argv) {
   }
   Grid_finalize();
 }  // main
+
+
+/*
+Input file example
+
+
+JSON
+
+{
+    "WilsonFlow":{
+	"steps": 200,
+	"step_size": 0.01,
+	"meas_interval": 50,
+  "maxTau": 2.0
+    },
+    "Configurations":{
+	"conf_prefix": "ckpoint_lat",
+	"rng_prefix": "ckpoint_rng",
+	"StartConfiguration": 3000,
+	"EndConfiguration": 3000,
+	"Skip": 5
+    }
+}
+
+
+*/