Merge pull request #412 from giltirn/patch/adaptive-wflow

Patch/adaptive wflow

Grid/lattice/Lattice_reduction.h

@@ -91,10 +91,7 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites)
   for(int i=0;i<nthread;i++){
     ssum = ssum+sumarray[i];
   } 
-  
-  typedef typename vobj::scalar_object ssobj;
-  ssobj ret = ssum;
-  return ret;
+  return ssum;
 }
 /*
 Threaded max, don't use for now

Grid/qcd/hmc/HMC.h

@@ -53,6 +53,7 @@ struct HMCparameters: Serializable {
                                   Integer, Trajectories, /* @brief Number of sweeps in this run */
                                   bool, MetropolisTest,
                                   Integer, NoMetropolisUntil,
+				  bool, PerformRandomShift, /* @brief Randomly shift the gauge configuration at the start of a trajectory */
                                   std::string, StartingType,
                                   IntegratorParameters, MD)
 
@@ -63,6 +64,7 @@ struct HMCparameters: Serializable {
     StartTrajectory   = 0;
     Trajectories      = 10;
     StartingType      = "HotStart";
+    PerformRandomShift = true;
     /////////////////////////////////
   }
 
@@ -83,6 +85,7 @@ struct HMCparameters: Serializable {
     std::cout << GridLogMessage << "[HMC parameters] Start trajectory        : " << StartTrajectory << "\n";
     std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
     std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs    : " << NoMetropolisUntil << "\n";
+    std::cout << GridLogMessage << "[HMC parameters] Doing random shift      : " << std::boolalpha << PerformRandomShift << "\n";
     std::cout << GridLogMessage << "[HMC parameters] Starting type           : " << StartingType << "\n";
     MD.print_parameters();
   }
@@ -95,6 +98,7 @@ private:
   const HMCparameters Params;
 
   typedef typename IntegratorType::Field Field;
+  typedef typename IntegratorType::FieldImplementation FieldImplementation;
   typedef std::vector< HmcObservable<Field> * > ObsListType;
 
   //pass these from the resource manager
@@ -138,26 +142,37 @@ private:
 
     GridBase *Grid = U.Grid();
 
-    //////////////////////////////////////////////////////////////////////////////////////////////////////
-    // Mainly for DDHMC perform a random translation of U modulo volume
-    //////////////////////////////////////////////////////////////////////////////////////////////////////
-    std::cout << GridLogMessage << "--------------------------------------------------\n";
-    std::cout << GridLogMessage << "Random shifting gauge field by [";
-    for(int d=0;d<Grid->Nd();d++) {
+    if(Params.PerformRandomShift){
+      //////////////////////////////////////////////////////////////////////////////////////////////////////
+      // Mainly for DDHMC perform a random translation of U modulo volume
+      //////////////////////////////////////////////////////////////////////////////////////////////////////
+      std::cout << GridLogMessage << "--------------------------------------------------\n";
+      std::cout << GridLogMessage << "Random shifting gauge field by [";
 
-      int L = Grid->GlobalDimensions()[d];
+      std::vector<typename FieldImplementation::GaugeLinkField> Umu(Grid->Nd(), U.Grid());
+      for(int mu=0;mu<Grid->Nd();mu++) Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
 
-      RealD rn_uniform;  random(sRNG, rn_uniform);
+      for(int d=0;d<Grid->Nd();d++) {
 
-      int shift = (int) (rn_uniform*L);
+	int L = Grid->GlobalDimensions()[d];
 
-      std::cout << shift;
-      if(d<Grid->Nd()-1) std::cout <<",";
-      else               std::cout <<"]\n";
+	RealD rn_uniform;  random(sRNG, rn_uniform);
+
+	int shift = (int) (rn_uniform*L);
+
+	std::cout << shift;
+	if(d<Grid->Nd()-1) std::cout <<",";
+	else               std::cout <<"]\n";
       
-      U = Cshift(U,d,shift);
+	//shift all fields together in a way that respects the gauge BCs
+	for(int mu=0; mu < Grid->Nd(); mu++)
+	  Umu[mu] = FieldImplementation::CshiftLink(Umu[mu],d,shift);
+      }
+
+      for(int mu=0;mu<Grid->Nd();mu++) PokeIndex<LorentzIndex>(U,Umu[mu],mu);
+      
+      std::cout << GridLogMessage << "--------------------------------------------------\n";
     }
-    std::cout << GridLogMessage << "--------------------------------------------------\n";
 
     TheIntegrator.reset_timer();
     

Grid/qcd/hmc/integrators/Integrator.h

@@ -63,10 +63,10 @@ public:
 };
 
 /*! @brief Class for Molecular Dynamics management */
-template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
+template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
 protected:
-
+  typedef FieldImplementation_ FieldImplementation;
   typedef typename FieldImplementation::Field MomentaField;  //for readability
   typedef typename FieldImplementation::Field Field;
 

Grid/qcd/hmc/integrators/Integrator_algorithm.h

@@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid);
  *  P 1/2                            P 1/2
  */
 
-template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
-class LeapFrog : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> 
+template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
+class LeapFrog : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy> 
 {
 public:
+  typedef FieldImplementation_ FieldImplementation;
   typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy> Algorithm;
   INHERIT_FIELD_TYPES(FieldImplementation);
 
@@ -135,13 +136,14 @@ public:
   }
 };
 
-template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
-class MinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> 
+template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
+class MinimumNorm2 : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy> 
 {
 private:
   const RealD lambda = 0.1931833275037836;
 
 public:
+  typedef FieldImplementation_ FieldImplementation;
   INHERIT_FIELD_TYPES(FieldImplementation);
 
   MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
@@ -192,8 +194,8 @@ public:
   }
 };
 
-template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
-class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> 
+template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
+class ForceGradient : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy> 
 {
 private:
   const RealD lambda = 1.0 / 6.0;
@@ -202,6 +204,7 @@ private:
   const RealD theta = 0.0;
 
 public:
+  typedef FieldImplementation_ FieldImplementation;
   INHERIT_FIELD_TYPES(FieldImplementation);
 
   // Looks like dH scales as dt^4. tested wilson/wilson 2 level.

Grid/qcd/observables/topological_charge.h

@@ -31,15 +31,16 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
+
 struct TopologySmearingParameters : Serializable {
   GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters,
-				  int, steps,
-				  float, step_size,
 				  int, meas_interval,
-				  float, maxTau);
+				  float, init_step_size,
+				  float, maxTau,
+				  float, tolerance);
 
-  TopologySmearingParameters(int s = 0, float ss = 0.0f, int mi = 0, float mT = 0.0f):
-    steps(s), step_size(ss), meas_interval(mi), maxTau(mT){}
+ TopologySmearingParameters(float ss = 0.0f, int mi = 0, float mT = 0.0f, float tol = 1e-4):
+  init_step_size(ss), meas_interval(mi), maxTau(mT), tolerance(tol){}
 
   template < class ReaderClass >
   TopologySmearingParameters(Reader<ReaderClass>& Reader){
@@ -97,8 +98,8 @@ public:
         
       if (Pars.do_smearing){
 	// using wilson flow by default here
-	WilsonFlow<PeriodicGimplR> WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval);
-	WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau);
+	WilsonFlowAdaptive<PeriodicGimplR> WF(Pars.Smearing.init_step_size, Pars.Smearing.maxTau, Pars.Smearing.tolerance, Pars.Smearing.meas_interval);
+	WF.smear(Usmear, U);
 	Real T0   = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear);
 	std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
 		  << "T0                : [ " << traj << " ] "<< T0 << std::endl;

Grid/qcd/smearing/WilsonFlow.h

@@ -33,27 +33,25 @@ directory
 NAMESPACE_BEGIN(Grid);
 
 template <class Gimpl>
-class WilsonFlow: public Smear<Gimpl>{
+class WilsonFlowBase: public Smear<Gimpl>{
 public:
   //Store generic measurements to take during smearing process using std::function
   typedef std::function<void(int, RealD, const typename Gimpl::GaugeField &)> FunctionType;  //int: step,  RealD: flow time,  GaugeField : the gauge field
-  
-private:
-  unsigned int Nstep;
-  RealD epsilon; //for regular smearing this is the time step, for adaptive it is the initial time step
- 
+
+protected:
   std::vector< std::pair<int, FunctionType> > functions; //The int maps to the measurement frequency
 
   mutable WilsonGaugeAction<Gimpl> SG;
-
-  //Evolve the gauge field by 1 step and update tau
-  void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const;
-  //Evolve the gauge field by 1 step and update tau and the current time step eps
-  void evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps, RealD maxTau) const;
-
+   
 public:
   INHERIT_GIMPL_TYPES(Gimpl)
 
+  explicit WilsonFlowBase(unsigned int meas_interval =1):
+    SG(WilsonGaugeAction<Gimpl>(3.0)) {
+    // WilsonGaugeAction with beta 3.0
+    setDefaultMeasurements(meas_interval);
+  }
+    
   void resetActions(){ functions.clear(); }
 
   void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); }
@@ -64,34 +62,11 @@ public:
   //and output to stdout
   void setDefaultMeasurements(int topq_meas_interval = 1);
 
-  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
-  Nstep(Nstep),
-    epsilon(epsilon),
-    SG(WilsonGaugeAction<Gimpl>(3.0)) {
-    // WilsonGaugeAction with beta 3.0
-    assert(epsilon > 0.0);
-    LogMessage();
-    setDefaultMeasurements(interval);
-  }
-
-  void LogMessage() {
-    std::cout << GridLogMessage
-	      << "[WilsonFlow] Nstep   : " << Nstep << std::endl;
-    std::cout << GridLogMessage
-	      << "[WilsonFlow] epsilon : " << epsilon << std::endl;
-    std::cout << GridLogMessage
-	      << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
-  }
-
-  virtual void smear(GaugeField&, const GaugeField&) const;
-
-  virtual void derivative(GaugeField&, const GaugeField&, const GaugeField&) const {
+  void derivative(GaugeField&, const GaugeField&, const GaugeField&) const override{
     assert(0);
     // undefined for WilsonFlow
   }
 
-  void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau) const;
-
   //Compute t^2 <E(t)> for time t from the plaquette
   static RealD energyDensityPlaquette(const RealD t, const GaugeField& U);
 
@@ -115,82 +90,63 @@ public:
   std::vector<RealD> flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1);
 };
 
+//Basic iterative Wilson flow
+template <class Gimpl>
+class WilsonFlow: public WilsonFlowBase<Gimpl>{
+private:
+  int Nstep; //number of steps
+  RealD epsilon;  //step size
+
+  //Evolve the gauge field by 1 step of size eps and update tau
+  void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const;
+
+public:
+  INHERIT_GIMPL_TYPES(Gimpl)
+
+  //Integrate the Wilson flow for Nstep steps of size epsilon
+  WilsonFlow(const RealD epsilon, const int Nstep, unsigned int meas_interval = 1): WilsonFlowBase<Gimpl>(meas_interval), Nstep(Nstep), epsilon(epsilon){}
+
+  void smear(GaugeField& out, const GaugeField& in) const override;
+};
+
+//Wilson flow with adaptive step size
+template <class Gimpl>
+class WilsonFlowAdaptive: public WilsonFlowBase<Gimpl>{
+private:
+  RealD init_epsilon; //initial step size
+  RealD maxTau; //integrate to t=maxTau
+  RealD tolerance; //integration error tolerance
+
+  //Evolve the gauge field by 1 step and update tau and the current time step eps
+  //
+  //If the step size eps is too large that a significant integration error results,
+  //the gauge field (U) and tau will not be updated and the function will return 0; eps will be adjusted to a smaller
+  //value for the next iteration.
+  //
+  //For a successful integration step the function will return 1
+  int evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps) const;
+
+public:
+  INHERIT_GIMPL_TYPES(Gimpl)
+
+  WilsonFlowAdaptive(const RealD init_epsilon, const RealD maxTau, const RealD tolerance, unsigned int meas_interval = 1): 
+  WilsonFlowBase<Gimpl>(meas_interval), init_epsilon(init_epsilon), maxTau(maxTau), tolerance(tolerance){}
+
+  void smear(GaugeField& out, const GaugeField& in) const override;
+};
 
 ////////////////////////////////////////////////////////////////////////////////
 // Implementations
 ////////////////////////////////////////////////////////////////////////////////
 template <class Gimpl>
-void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{
-  GaugeField Z(U.Grid());
-  GaugeField tmp(U.Grid());
-  SG.deriv(U, 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
-  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
-  tau += epsilon;
-}
-
-template <class Gimpl>
-void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps, RealD maxTau) const{
-  if (maxTau - tau < eps){
-    eps = maxTau-tau;
-  }
-  //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*eps);    // 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*eps);    // 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*eps);    // V(t+e) = exp(ep*Z)*W2
-
-  // Ramos 
-  Gimpl::update_field(Zprime, Uprime, -2.0*eps); // 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
-    
-  tau += eps;
-  //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
-    
-  eps = eps*0.95*std::pow(1e-4/diff,1./3.);
-  //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
-
-}
-
-
-template <class Gimpl>
-RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
+RealD WilsonFlowBase<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
   static WilsonGaugeAction<Gimpl> SG(3.0);
   return 2.0 * t * t * SG.S(U)/U.Grid()->gSites();
 }
 
 //Compute t^2 <E(t)> for time from the 1x1 cloverleaf form
 template <class Gimpl>
-RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){
+RealD WilsonFlowBase<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){
   typedef typename Gimpl::GaugeLinkField GaugeMat;
   typedef typename Gimpl::GaugeField GaugeLorentz;
 
@@ -215,7 +171,7 @@ RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField
 
 
 template <class Gimpl>
-std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){
+std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){
   std::vector<RealD> out;
   resetActions();
   addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ 
@@ -227,13 +183,13 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeFie
 }
 
 template <class Gimpl>
-std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){
+std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){
   GaugeField V(U);
   return flowMeasureEnergyDensityPlaquette(V,U, measure_interval);
 }
 
 template <class Gimpl>
-std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){
+std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){
   std::vector<RealD> out;
   resetActions();
   addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ 
@@ -245,16 +201,52 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeFi
 }
 
 template <class Gimpl>
-std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){
+std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){
   GaugeField V(U);
   return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval);
 }
 
+template <class Gimpl>
+void WilsonFlowBase<Gimpl>::setDefaultMeasurements(int topq_meas_interval){
+  addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){
+      std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "  << step << "  " << t << "  " << energyDensityPlaquette(t,U) << std::endl;
+    });
+  addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
+      std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "  << step << "  " << WilsonLoops<Gimpl>::TopologicalCharge(U) << std::endl;
+    });
+}
 
 
-//#define WF_TIMING 
+
+template <class Gimpl>
+void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{
+  GaugeField Z(U.Grid());
+  GaugeField tmp(U.Grid());
+  this->SG.deriv(U, 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;
+  this->SG.deriv(U, tmp); Z += tmp;                 // -17/32*Z0 +Z1
+  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;
+  this->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
+  tau += epsilon;
+}
+
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] Nstep   : " << Nstep << std::endl;
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] epsilon : " << epsilon << std::endl;
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
+
   out = in;
   RealD taus = 0.;
   for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement
@@ -266,37 +258,93 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
     std::cout << "Time to evolve " << diff.count() << " s\n";
 #endif
     //Perform measurements
-    for(auto const &meas : functions)
+    for(auto const &meas : this->functions)
       if( step % meas.first == 0 ) meas.second(step,taus,out);
   }
 }
 
+
+
 template <class Gimpl>
-void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau) const{
-  out = in;
-  RealD taus = 0.;
-  RealD eps = epsilon;
-  unsigned int step = 0;
-  do{
-    step++;
-    //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
-    evolve_step_adaptive(out, taus, eps, maxTau);
-    //Perform measurements
-    for(auto const &meas : functions)
-      if( step % meas.first == 0 ) meas.second(step,taus,out);
-  } while (taus < maxTau);
+int WilsonFlowAdaptive<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps) const{
+  if (maxTau - tau < eps){
+    eps = maxTau-tau;
+  }
+  //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
+  GaugeField Z(U.Grid());
+  GaugeField Zprime(U.Grid());
+  GaugeField tmp(U.Grid()), Uprime(U.Grid()), Usave(U.Grid());
+  Uprime = U;
+  Usave = U;
+
+  this->SG.deriv(U, Z);
+  Zprime = -Z;
+  Z *= 0.25;                                  // Z0 = 1/4 * F(U)
+  Gimpl::update_field(Z, U, -2.0*eps);    // U = W1 = exp(ep*Z0)*W0
+
+  Z *= -17.0/8.0;
+  this->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*eps);    // U_= W2 = exp(ep*Z)*W1
+    
+
+  Z *= -4.0/3.0;
+  this->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*eps);    // V(t+e) = exp(ep*Z)*W2
+
+  // Ramos arXiv:1301.4388
+  Gimpl::update_field(Zprime, Uprime, -2.0*eps); // V'(t+e) = exp(ep*Z')*W0
+
+  // Compute distance using Ramos' definition
+  GaugeField diffU = U - Uprime;
+  RealD max_dist = 0;
+
+  for(int mu=0;mu<Nd;mu++){
+    typename Gimpl::GaugeLinkField diffU_mu = PeekIndex<LorentzIndex>(diffU, mu);
+    RealD dist_mu = sqrt( maxLocalNorm2(diffU_mu) ) /Nc/Nc; //maximize over sites
+    max_dist = std::max(max_dist, dist_mu); //maximize over mu
+  }
+  
+  int ret;
+  if(max_dist < tolerance) {
+    tau += eps;
+    ret = 1;
+  } else {
+    U = Usave;
+    ret = 0;
+  }
+  eps = eps*0.95*std::pow(tolerance/max_dist,1./3.);
+  std::cout << GridLogMessage << "Adaptive smearing : Distance: "<< max_dist <<" Step successful: " << ret << " New epsilon: " << eps << std::endl; 
+
+  return ret;
 }
 
 template <class Gimpl>
-void WilsonFlow<Gimpl>::setDefaultMeasurements(int topq_meas_interval){
-  addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){
-      std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "  << step << "  " << t << "  " << energyDensityPlaquette(t,U) << std::endl;
-    });
-  addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
-      std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "  << step << "  " << WilsonLoops<Gimpl>::TopologicalCharge(U) << std::endl;
-    });
+void WilsonFlowAdaptive<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] initial epsilon : " << init_epsilon << std::endl;
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] full trajectory : " << maxTau << std::endl;
+  std::cout << GridLogMessage
+	    << "[WilsonFlow] tolerance   : " << tolerance << std::endl;
+  out = in;
+  RealD taus = 0.;
+  RealD eps = init_epsilon;
+  unsigned int step = 0;
+  do{
+    int step_success = evolve_step_adaptive(out, taus, eps); 
+    step += step_success; //step will not be incremented if the integration step fails
+
+    //Perform measurements
+    if(step_success)
+      for(auto const &meas : this->functions)
+	if( step % meas.first == 0 ) meas.second(step,taus,out);
+  } while (taus < maxTau);
 }
 
 
+
 NAMESPACE_END(Grid);
 

Grid/qcd/utils/CovariantCshift.h

@@ -227,26 +227,38 @@ namespace ConjugateBC {
   //shift = -1
   //Out(x) = U_\mu(x-mu)  | x_\mu != 0
   //       = U*_\mu(L-1)  | x_\mu == 0
+  //shift = 2
+  //Out(x) = U_\mu(x+2\hat\mu)  | x_\mu < L-2
+  //       = U*_\mu(1)  | x_\mu == L-1
+  //       = U*_\mu(0)  | x_\mu == L-2
+  //shift = -2
+  //Out(x) = U_\mu(x-2mu)  | x_\mu > 1
+  //       = U*_\mu(L-2)  | x_\mu == 0
+  //       = U*_\mu(L-1)  | x_\mu == 1
+  //etc
   template<class gauge> Lattice<gauge>
   CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
   {
     GridBase *grid = Link.Grid();
-    int Lmu = grid->GlobalDimensions()[mu] - 1;
+    int Lmu = grid->GlobalDimensions()[mu];
+    assert(abs(shift) < Lmu && "Invalid shift value");
 
     Lattice<iScalar<vInteger>> coor(grid);
     LatticeCoordinate(coor, mu);
 
     Lattice<gauge> tmp(grid);
-    if(shift == 1){
-      tmp = Cshift(Link, mu, 1);
-      tmp = where(coor == Lmu, conjugate(tmp), tmp);
+    if(shift > 0){
+      tmp = Cshift(Link, mu, shift);
+      tmp = where(coor >= Lmu-shift, conjugate(tmp), tmp);
       return tmp;
-    }else if(shift == -1){
+    }else if(shift < 0){
       tmp = Link;
-      tmp = where(coor == Lmu, conjugate(tmp), tmp);
-      return Cshift(tmp, mu, -1);
-    }else assert(0 && "Invalid shift value");
-    return tmp; //shuts up the compiler fussing about the return type
+      tmp = where(coor >= Lmu+shift, conjugate(tmp), tmp);
+      return Cshift(tmp, mu, shift);
+    }
+    
+    //shift == 0
+    return Link;
   }
 
 }

Grid/qcd/utils/GaugeFix.h

@@ -72,12 +72,12 @@ public:
 
   //Fix the gauge field Umu
   //0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
     GridBase *grid = Umu.Grid();
     GaugeMat xform(grid);
     SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge);
   }
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
   //Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform
 
     GridBase *grid = Umu.Grid();

Grid/simd/Grid_sse4.h

@@ -35,7 +35,7 @@ Author: neo <cossu@post.kek.jp>
 */
 // Time-stamp: <2015-06-16 23:27:54 neo>
 //----------------------------------------------------------------------
-
+#include <immintrin.h>
 #include <pmmintrin.h>
 
 NAMESPACE_BEGIN(Grid);

tests/Test_gfield_shift.cc

@@ -0,0 +1,183 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_gfield_shift.cc
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@bnl.gov>
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+
+//Test the shifting of the gauge field that respects the boundary conditions
+#include <Grid/Grid.h>
+
+using namespace Grid;
+ ;
+
+typedef ConjugateGimplR Gimpl; //can choose periodic / charge conjugate directions at wil
+typedef Gimpl::GaugeField GaugeField;
+typedef Gimpl::GaugeLinkField GaugeLinkField;
+typedef Gimpl::SiteGaugeField SiteGaugeField;
+typedef Gimpl::SiteGaugeLink SiteGaugeLink;
+
+GaugeField CshiftGaugeField(const GaugeField &U, const int dir, const int shift){
+  GridBase *Grid = U.Grid();
+
+  GaugeField out(Grid);
+  GaugeLinkField Umu(Grid);
+  for(int mu=0;mu<Grid->Nd();mu++){
+    Umu = PeekIndex<LorentzIndex>(U, mu);
+    Umu = Gimpl::CshiftLink(Umu,dir,shift);
+    PokeIndex<LorentzIndex>(out,Umu,mu);
+  }
+  return out;
+}
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+
+  std::vector<int> conj_dirs = {1,1,0,0}; 
+  Gimpl::setDirections(conj_dirs);
+
+  GridCartesian        Fine(latt_size,simd_layout,mpi_layout);
+
+  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+
+  GaugeField U(&Fine);
+  GaugeField ShiftU(&Fine);
+
+  GaugeLinkField link_field(&Fine), link_field_2(&Fine);
+
+  //Like Test_cshift we put the lex coordinate index on each site but make it imaginary
+  //so we can tell when it was complex conjugated
+  LatticeComplex lex(&Fine);
+  lex=Zero();
+  U = Zero();
+  {
+    LatticeComplex coor(&Fine);
+    Integer stride =1;
+    for(int d=0;d<4;d++){
+      LatticeCoordinate(coor,d);
+      lex = lex + coor*stride;
+      stride=stride*latt_size[d];
+    }
+    PokeIndex<ColourIndex>(link_field, lex, 0,0); //place on 0,0 element of link
+
+    for(int mu=0;mu<Nd;mu++){
+      link_field_2 = link_field + mu*stride; //add in lex-mapping of mu
+      link_field_2 = ComplexD(0,1) * link_field_2; //make imaginary
+      PokeIndex<LorentzIndex>(U, link_field_2, mu);
+    }
+  }
+
+  std::stringstream ss;
+  ss<<"error";
+  for(int d=0;d<Fine._ndimension;d++){
+    ss<<"."<<Fine._processor_coor[d];
+  }
+  ss<<"_wr_"<<Fine._processor;
+  std::string fname(ss.str());
+  std::ofstream ferr(fname);
+  
+  Integer vol4d = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+
+  bool fail = false;
+  typename SiteGaugeField::scalar_object um;
+  TComplex cm;
+  
+  for(int dir=0;dir<4;dir++){
+    for(int shift=-latt_size[dir]+1;shift<latt_size[dir];shift++){
+      if ( Fine.IsBoss() )
+	std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction "<<dir
+		 << " dir is conj ? " << conj_dirs[dir] << std::endl;
+
+      ShiftU = CshiftGaugeField(U,dir,shift);
+
+      Coordinate coor(4);
+      
+      for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
+      for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
+      for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
+      for(coor[0]=0;coor[0]<latt_size[0];coor[0]++){  
+	peekSite(um,ShiftU,coor);
+
+	Coordinate scoor(coor);
+	scoor[dir] = (scoor[dir]+shift + latt_size[dir])%latt_size[dir];
+	
+	Integer slex = scoor[0]
+	  + latt_size[0]*scoor[1]
+	  + latt_size[0]*latt_size[1]*scoor[2]
+	  + latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
+	
+	for(int mu = 0 ; mu < 4; mu++){
+	  Integer slex_mu = slex + vol4d*mu;
+	  Complex scm(0,slex_mu); //imaginary
+	  if(
+	     ( shift > 0 && coor[dir] >= latt_size[dir]-shift && conj_dirs[dir] ) 
+	     ||
+	     ( shift < 0 && coor[dir] <= -shift-1 && conj_dirs[dir] )
+	     )
+	    scm = conjugate(scm); //CC if pulled over boundary
+
+	  cm = um(mu)()(0,0);
+
+	  RealD nrm = abs(scm-cm()()());
+	  //std::cout << cm << " " << scm << std::endl;
+
+	  Coordinate peer(4);
+	  Complex tmp  =cm;
+	  Integer index=real(tmp);
+
+	  Integer cm_mu = index / vol4d;
+	  index = index % vol4d;
+	  Lexicographic::CoorFromIndex(peer,index,latt_size);
+
+	  if (nrm > 0){
+	    ferr<<"FAIL mu " << mu << " shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<"  "<<nrm<<std::endl;
+	    ferr<<"Got mu "<< cm_mu << " site " <<index<<" : " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
+
+	    index=real(scm);
+	    Integer scm_mu = index / vol4d;
+	    index = index % vol4d;
+	    Lexicographic::CoorFromIndex(peer,index,latt_size);
+	    ferr<<"Expect mu " << scm_mu << " site " <<index<<": " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
+	    fail = true;
+	  }
+	}
+
+	}}}}
+    }
+  }
+
+  if(fail) std::cout << "Test FAILED : see " << fname << " for more details" << std::endl;
+  else std::cout << "Test Passed" << std::endl;
+
+  Grid_finalize();
+}

tests/smearing/Test_WilsonFlow_adaptive.cc

@@ -0,0 +1,153 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/hmc/Test_WilsonFlow_adaptive.cc
+
+Copyright (C) 2017
+
+Author: Christopher Kelly <ckelly@bnl.gov>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+//Linearly interpolate between two nearest times
+RealD interpolate(const RealD t_int, const std::vector<std::pair<RealD,RealD> > &data){
+  RealD tdiff1=1e32; int t1_idx=-1;
+  RealD tdiff2=1e32; int t2_idx=-1;
+
+  for(int i=0;i<data.size();i++){
+    RealD diff = fabs(data[i].first-t_int);
+    //std::cout << "targ " << t_int << " cur " << data[i].first << " diff " << diff << " best diff1 " << tdiff1 << " diff2 " << tdiff2 << std::endl;
+
+    if(diff < tdiff1){ 
+      if(tdiff1 < tdiff2){ //swap out tdiff2
+	tdiff2 = tdiff1; t2_idx = t1_idx;
+      }
+      tdiff1 = diff; t1_idx = i; 
+    }
+    else if(diff < tdiff2){ tdiff2 = diff; t2_idx = i; }
+  }
+  assert(t1_idx != -1 && t2_idx != -1);
+  
+  RealD t2 = data[t2_idx].first,  v2 = data[t2_idx].second;
+  RealD t1 = data[t1_idx].first,  v1 = data[t1_idx].second;
+  
+  //v = a + bt
+  //v2-v1 = b(t2-t1)
+  RealD b = (v2-v1)/(t2-t1);
+  RealD a = v1 - b*t1;
+  RealD vout = a + b*t_int;
+
+  //std::cout << "Interpolate to " << t_int << " two closest points " << t1  << " " << t2 
+  //<< " with values " << v1 << " "<< v2 << " : got " << vout <<  std::endl;
+  return vout;
+}
+
+
+int main(int argc, char **argv) {
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
+
+  std::vector<int> seeds({1, 2, 3, 4, 5});
+  GridSerialRNG sRNG;
+  GridParallelRNG pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+
+  LatticeGaugeField U(&Grid);
+  SU<Nc>::HotConfiguration(pRNG, U);
+
+  int Nstep = 300;
+  RealD epsilon = 0.01;
+  RealD maxTau = Nstep*epsilon;
+  RealD tolerance = 1e-4;
+
+  for(int i=1;i<argc;i++){
+    std::string sarg(argv[i]);
+    if(sarg == "--tolerance"){
+      std::stringstream ss; ss << argv[i+1]; ss >> tolerance;
+    }
+  }
+  std::cout << "Adaptive smear tolerance " << tolerance << std::endl;
+
+  //Setup iterative Wilson flow
+  WilsonFlow<PeriodicGimplD> wflow(epsilon,Nstep);
+  wflow.resetActions();    
+
+  std::vector<std::pair<RealD, RealD> > meas_orig;
+
+  wflow.addMeasurement(1, [&wflow,&meas_orig](int step, RealD t, const LatticeGaugeField &U){ 
+      std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl;
+      meas_orig.push_back( {t, wflow.energyDensityCloverleaf(t,U)} );
+    });
+
+  //Setup adaptive Wilson flow
+  WilsonFlowAdaptive<PeriodicGimplD> wflow_ad(epsilon,maxTau,tolerance);
+  wflow_ad.resetActions();    
+
+  std::vector<std::pair<RealD, RealD> > meas_adaptive;
+
+  wflow_ad.addMeasurement(1, [&wflow_ad,&meas_adaptive](int step, RealD t, const LatticeGaugeField &U){ 
+      std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl;
+      meas_adaptive.push_back( {t, wflow_ad.energyDensityCloverleaf(t,U)} );
+    });
+
+  //Run
+  LatticeGaugeFieldD Vtmp(U.Grid());
+  wflow.smear(Vtmp, U); //basic smear
+
+  Vtmp = Zero();
+  wflow_ad.smear(Vtmp, U);
+
+  //Output values for plotting
+  {
+    std::ofstream out("wflow_t2E_orig.dat");
+    out.precision(16);
+    for(auto const &e: meas_orig){
+      out << e.first << " " << e.second << std::endl;
+    }
+  }
+  {
+    std::ofstream out("wflow_t2E_adaptive.dat");
+    out.precision(16);
+    for(auto const &e: meas_adaptive){
+      out << e.first << " " << e.second << std::endl;
+    }
+  }
+  
+  //Compare at times available with adaptive smearing
+  for(int i=0;i<meas_adaptive.size();i++){
+    RealD t = meas_adaptive[i].first;
+    RealD v_adaptive = meas_adaptive[i].second;
+    RealD v_orig = interpolate(t,  meas_orig); //should be very precise due to fine timestep
+    std::cout << t << " orig: " << v_orig << " adaptive: " << v_adaptive << " reldiff: " << (v_adaptive-v_orig)/v_orig << std::endl;
+  }
+  
+  std::cout << GridLogMessage << "Done" << std::endl;
+  Grid_finalize();
+}