HMC running with new formulation

Grid/lattice/Lattice.h

@@ -46,3 +46,4 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_unary.h>
 #include <Grid/lattice/Lattice_transfer.h>
 #include <Grid/lattice/Lattice_basis.h>
+#include <Grid/lattice/Lattice_crc.h>

Grid/lattice/Lattice_crc.h

@@ -0,0 +1,55 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/lattice/Lattice_crc.h
+
+    Copyright (C) 2021
+
+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 */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+template<class vobj> void DumpSliceNorm(std::string s,Lattice<vobj> &f,int mu=-1)
+{
+  auto ff = localNorm2(f);
+  if ( mu==-1 ) mu = f.Grid()->Nd()-1;
+  typedef typename vobj::tensor_reduced normtype;
+  typedef typename normtype::scalar_object scalar;
+  std::vector<scalar> sff;
+  sliceSum(ff,sff,mu);
+  for(int t=0;t<sff.size();t++){
+    std::cout << s<<" "<<t<<" "<<sff[t]<<std::endl;
+  }
+}
+
+template<class vobj> uint32_t crc(Lattice<vobj> & buf)
+{
+  autoView( buf_v , buf, CpuRead);
+  return ::crc32(0L,(unsigned char *)&buf_v[0],(size_t)sizeof(vobj)*buf.oSites());
+}
+
+#define CRC(U) std::cout << "FingerPrint "<<__FILE__ <<" "<< __LINE__ <<" "<< #U <<" "<<crc(U)<<std::endl;
+
+NAMESPACE_END(Grid);
+
+

Grid/log/Log.cc

@@ -69,6 +69,7 @@ GridLogger GridLogDebug  (1, "Debug", GridLogColours, "PURPLE");
 GridLogger GridLogPerformance(1, "Performance", GridLogColours, "GREEN");
 GridLogger GridLogIterative  (1, "Iterative", GridLogColours, "BLUE");
 GridLogger GridLogIntegrator (1, "Integrator", GridLogColours, "BLUE");
+GridLogger GridLogHMC (1, "HMC", GridLogColours, "BLUE");
 
 void GridLogConfigure(std::vector<std::string> &logstreams) {
   GridLogError.Active(0);
@@ -79,6 +80,7 @@ void GridLogConfigure(std::vector<std::string> &logstreams) {
   GridLogPerformance.Active(0);
   GridLogIntegrator.Active(1);
   GridLogColours.Active(0);
+  GridLogHMC.Active(1);
 
   for (int i = 0; i < logstreams.size(); i++) {
     if (logstreams[i] == std::string("Error"))       GridLogError.Active(1);
@@ -87,7 +89,8 @@ void GridLogConfigure(std::vector<std::string> &logstreams) {
     if (logstreams[i] == std::string("Iterative"))   GridLogIterative.Active(1);
     if (logstreams[i] == std::string("Debug"))       GridLogDebug.Active(1);
     if (logstreams[i] == std::string("Performance")) GridLogPerformance.Active(1);
-    if (logstreams[i] == std::string("Integrator"))  GridLogIntegrator.Active(1);
+    if (logstreams[i] == std::string("NoIntegrator"))  GridLogIntegrator.Active(0);
+    if (logstreams[i] == std::string("NoHMC"))         GridLogHMC.Active(0);
     if (logstreams[i] == std::string("Colours"))     GridLogColours.Active(1);
   }
 }

Grid/log/Log.h

@@ -182,6 +182,7 @@ extern GridLogger GridLogDebug  ;
 extern GridLogger GridLogPerformance;
 extern GridLogger GridLogIterative  ;
 extern GridLogger GridLogIntegrator  ;
+extern GridLogger GridLogHMC;
 extern Colours    GridLogColours;
 
 std::string demangle(const char* name) ;

Grid/qcd/action/ActionBase.h

@@ -40,6 +40,29 @@ class Action
 
 public:
   bool is_smeared = false;
+  RealD deriv_norm_sum;
+  RealD deriv_max_sum;
+  int   deriv_num;
+  RealD deriv_us;
+  RealD S_us;
+  RealD refresh_us;
+  void  reset_timer(void)        {
+    deriv_us = S_us = refresh_us = 0.0;
+    deriv_num=0;
+    deriv_norm_sum = deriv_max_sum=0.0;
+  }
+  void  deriv_log(RealD nrm, RealD max) { deriv_max_sum+=max; deriv_norm_sum+=nrm; deriv_num++;}
+  RealD deriv_max_average(void)         { return deriv_max_sum/deriv_num; };
+  RealD deriv_norm_average(void)        { return deriv_norm_sum/deriv_num; };
+  RealD deriv_timer(void)        { return deriv_us; };
+  RealD S_timer(void)            { return deriv_us; };
+  RealD refresh_timer(void)      { return deriv_us; };
+  void deriv_timer_start(void)   { deriv_us-=usecond(); }
+  void deriv_timer_stop(void)    { deriv_us+=usecond(); }
+  void refresh_timer_start(void) { refresh_us-=usecond(); }
+  void refresh_timer_stop(void)  { refresh_us+=usecond(); }
+  void S_timer_start(void)       { S_us-=usecond(); }
+  void S_timer_stop(void)        { S_us+=usecond(); }
   // Heatbath?
   virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
   virtual RealD S(const GaugeField& U) = 0;                             // evaluate the action

Grid/qcd/action/ActionCore.h

@@ -39,6 +39,7 @@ NAMESPACE_CHECK(ActionParams);
 
 #include <Grid/qcd/action/filters/MomentumFilter.h>
 #include <Grid/qcd/action/filters/DirichletFilter.h>
+#include <Grid/qcd/action/filters/DDHMCFilter.h>
 
 ////////////////////////////////////////////
 // Gauge Actions

Grid/qcd/hmc/GenericHMCrunner.h

@@ -129,18 +129,10 @@ public:
     Runner(S);
   }
 
-  //////////////////////////////////////////////////////////////////
+  //Use the checkpointer to initialize the RNGs and the gauge field, writing the resulting gauge field into U.
-
+  //This is called automatically by Run but may be useful elsewhere, e.g. for integrator tuning experiments
-private:
+  void initializeGaugeFieldAndRNGs(Field &U){
-  template <class SmearingPolicy>
+    if(!Resources.haveRNGs()) Resources.AddRNGs();
-  void Runner(SmearingPolicy &Smearing) {
-    auto UGrid = Resources.GetCartesian();
-    Resources.AddRNGs();
-    Field U(UGrid);
-
-    // Can move this outside?
-    typedef IntegratorType<SmearingPolicy> TheIntegrator;
-    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing);
 
     if (Parameters.StartingType == "HotStart") {
       // Hot start
@@ -167,6 +159,25 @@ private:
 	<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
       exit(1);
     }
+  }
+
+
+
+  //////////////////////////////////////////////////////////////////
+
+private:
+  template <class SmearingPolicy>
+  void Runner(SmearingPolicy &Smearing) {
+    auto UGrid = Resources.GetCartesian();
+    Field U(UGrid);
+
+    initializeGaugeFieldAndRNGs(U);
+
+    typedef IntegratorType<SmearingPolicy> TheIntegrator;
+    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing);
+
+    // Sets the momentum filter
+    MDynamics.setMomentumFilter(*(Resources.GetMomentumFilter()));
 
     Smearing.set_Field(U);
 

Grid/qcd/hmc/HMC.h

@@ -34,6 +34,7 @@ directory
 			    * @brief Classes for Hybrid Monte Carlo update
 			    *
 			    * @author Guido Cossu
+			    * @author Peter Boyle
 			    */
 			   //--------------------------------------------------------------------
 #pragma once
@@ -115,22 +116,17 @@ private:
 
     random(sRNG, rn_test);
 
-    std::cout << GridLogMessage
+    std::cout << GridLogHMC << "--------------------------------------------------\n";
-              << "--------------------------------------------------\n";
+    std::cout << GridLogHMC << "exp(-dH) = " << prob << "  Random = " << rn_test << "\n";
-    std::cout << GridLogMessage << "exp(-dH) = " << prob
+    std::cout << GridLogHMC << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
-              << "  Random = " << rn_test << "\n";
-    std::cout << GridLogMessage
-              << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
 
     if ((prob > 1.0) || (rn_test <= prob)) {  // accepted
-      std::cout << GridLogMessage << "Metropolis_test -- ACCEPTED\n";
+      std::cout << GridLogHMC << "Metropolis_test -- ACCEPTED\n";
-      std::cout << GridLogMessage
+      std::cout << GridLogHMC << "--------------------------------------------------\n";
-                << "--------------------------------------------------\n";
       return true;
     } else {  // rejected
-      std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n";
+      std::cout << GridLogHMC << "Metropolis_test -- REJECTED\n";
-      std::cout << GridLogMessage
+      std::cout << GridLogHMC << "--------------------------------------------------\n";
-                << "--------------------------------------------------\n";
       return false;
     }
   }
@@ -139,19 +135,68 @@ private:
   // Evolution
   /////////////////////////////////////////////////////////
   RealD evolve_hmc_step(Field &U) {
-    TheIntegrator.refresh(U, sRNG, pRNG);  // set U and initialize P and phi's
 
-    RealD H0 = TheIntegrator.S(U);  // initial state action
+    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++) {
+
+      int L = Grid->GlobalDimensions()[d];
+
+      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);
+    }
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+
+    TheIntegrator.reset_timer();
+    
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    // set U and initialize P and phi's
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << "Refresh momenta and pseudofermions";
+    TheIntegrator.refresh(U, sRNG, pRNG);  
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    // initial state action
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << "Compute initial action";
+    RealD H0 = TheIntegrator.S(U);  
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
 
     std::streamsize current_precision = std::cout.precision();
     std::cout.precision(15);
-    std::cout << GridLogMessage << "Total H before trajectory = " << H0 << "\n";
+    std::cout << GridLogHMC << "Total H before trajectory = " << H0 << "\n";
     std::cout.precision(current_precision);
 
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << " Molecular Dynamics evolution ";
     TheIntegrator.integrate(U);
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
 
-    RealD H1 = TheIntegrator.S(U);  // updated state action
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    // updated state action
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << "Compute final action";
+    RealD H1 = TheIntegrator.S(U);  
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
 
+
+    
     ///////////////////////////////////////////////////////////
     if(0){
       std::cout << "------------------------- Reversibility test" << std::endl;
@@ -163,17 +208,16 @@ private:
     }
     ///////////////////////////////////////////////////////////
 
-
     std::cout.precision(15);
-    std::cout << GridLogMessage << "Total H after trajectory  = " << H1
+
-	      << "  dH = " << H1 - H0 << "\n";
+    std::cout << GridLogHMC << "--------------------------------------------------\n";
+    std::cout << GridLogHMC << "Total H after trajectory  = " << H1 << "  dH = " << H1 - H0 << "\n";
+    std::cout << GridLogHMC << "--------------------------------------------------\n";
+
     std::cout.precision(current_precision);
     
     return (H1 - H0);
   }
-  
-
-  
 
 public:
   /////////////////////////////////////////
@@ -195,10 +239,13 @@ public:
 
     // Actual updates (evolve a copy Ucopy then copy back eventually)
     unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory;
+
     for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) {
-      std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
+
+      std::cout << GridLogHMC << "-- # Trajectory = " << traj << "\n";
+
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
-      	std::cout << GridLogMessage << "-- Thermalization" << std::endl;
+      	std::cout << GridLogHMC << "-- Thermalization" << std::endl;
       }
       
       double t0=usecond();
@@ -207,20 +254,19 @@ public:
       DeltaH = evolve_hmc_step(Ucopy);
       // Metropolis-Hastings test
       bool accept = true;
-      if (traj >= Params.StartTrajectory + Params.NoMetropolisUntil) {
+      if (Params.MetropolisTest && traj >= Params.StartTrajectory + Params.NoMetropolisUntil) {
         accept = metropolis_test(DeltaH);
       } else {
-      	std::cout << GridLogMessage << "Skipping Metropolis test" << std::endl;
+      	std::cout << GridLogHMC << "Skipping Metropolis test" << std::endl;
       }
 
       if (accept)
         Ucur = Ucopy; 
       
-     
-      
       double t1=usecond();
-      std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
+      std::cout << GridLogHMC << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 
+      TheIntegrator.print_timer();
 
       for (int obs = 0; obs < Observables.size(); obs++) {
       	std::cout << GridLogDebug << "Observables # " << obs << std::endl;
@@ -228,7 +274,7 @@ public:
       	std::cout << GridLogDebug << "Observables pointer " << Observables[obs] << std::endl;
         Observables[obs]->TrajectoryComplete(traj + 1, Ucur, sRNG, pRNG);
       }
-      std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::::" << std::endl;
+      std::cout << GridLogHMC << ":::::::::::::::::::::::::::::::::::::::::::" << std::endl;
     }
   }
 

Grid/qcd/hmc/HMCResourceManager.h

@@ -72,6 +72,8 @@ class HMCResourceManager {
   typedef HMCModuleBase< BaseHmcCheckpointer<ImplementationPolicy> > CheckpointerBaseModule;
   typedef HMCModuleBase< HmcObservable<typename ImplementationPolicy::Field> > ObservableBaseModule;
   typedef ActionModuleBase< Action<typename ImplementationPolicy::Field>, GridModule > ActionBaseModule;
+  typedef typename ImplementationPolicy::Field MomentaField;
+  typedef typename ImplementationPolicy::Field Field;  
 
   // Named storage for grid pairs (std + red-black)
   std::unordered_map<std::string, GridModule> Grids;
@@ -80,6 +82,9 @@ class HMCResourceManager {
   // SmearingModule<ImplementationPolicy> Smearing;
   std::unique_ptr<CheckpointerBaseModule> CP;
 
+  // Momentum filter
+  std::unique_ptr<MomentumFilterBase<typename ImplementationPolicy::Field> > Filter;
+  
   // A vector of HmcObservable modules
   std::vector<std::unique_ptr<ObservableBaseModule> > ObservablesList;
 
@@ -90,6 +95,7 @@ class HMCResourceManager {
 
   bool have_RNG;
   bool have_CheckPointer;
+  bool have_Filter;
 
   // NOTE: operator << is not overloaded for std::vector<string> 
   // so this function is necessary
@@ -101,7 +107,7 @@ class HMCResourceManager {
 
 
 public:
-  HMCResourceManager() : have_RNG(false), have_CheckPointer(false) {}
+  HMCResourceManager() : have_RNG(false), have_CheckPointer(false), have_Filter(false) {}
 
   template <class ReaderClass, class vector_type = vComplex >
   void initialize(ReaderClass &Read){
@@ -129,6 +135,7 @@ public:
     RNGModuleParameters RNGpar(Read);
     SetRNGSeeds(RNGpar);
 
+  
     // Observables
     auto &ObsFactory = HMC_ObservablesModuleFactory<observable_string, typename ImplementationPolicy::Field, ReaderClass>::getInstance(); 
     Read.push(observable_string);// here must check if existing...
@@ -208,6 +215,16 @@ public:
     AddGrid(s, Mod);
   }
 
+  void SetMomentumFilter( MomentumFilterBase<typename ImplementationPolicy::Field> * MomFilter) {
+    assert(have_Filter==false);
+    Filter = std::unique_ptr<MomentumFilterBase<typename ImplementationPolicy::Field> >(MomFilter);
+    have_Filter = true;
+  }
+  MomentumFilterBase<typename ImplementationPolicy::Field> *GetMomentumFilter(void) {
+    if ( !have_Filter)
+      SetMomentumFilter(new MomentumFilterNone<typename ImplementationPolicy::Field>());
+    return Filter.get();
+  }
 
   GridCartesian* GetCartesian(std::string s = "") {
     if (s.empty()) s = Grids.begin()->first;
@@ -226,6 +243,9 @@ public:
   //////////////////////////////////////////////////////
   // Random number generators
   //////////////////////////////////////////////////////
+  
+  //Return true if the RNG objects have been instantiated
+  bool haveRNGs() const{ return have_RNG; }
 
   void AddRNGs(std::string s = "") {
     // Couple the RNGs to the GridModule tagged by s

Grid/qcd/hmc/integrators/Integrator.h

@@ -66,6 +66,7 @@ public:
 template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
 protected:
+
   typedef typename FieldImplementation::Field MomentaField;  //for readability
   typedef typename FieldImplementation::Field Field;
 
@@ -118,36 +119,58 @@ protected:
     }
   } update_P_hireps{};
 
+ 
   void update_P(MomentaField& Mom, Field& U, int level, double ep) {
     // input U actually not used in the fundamental case
     // Fundamental updates, include smearing
 
     for (int a = 0; a < as[level].actions.size(); ++a) {
+
       double start_full = usecond();
       Field force(U.Grid());
       conformable(U.Grid(), Mom.Grid());
 
       Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
       double start_force = usecond();
+      as[level].actions.at(a)->deriv_timer_start();
       as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta
+      as[level].actions.at(a)->deriv_timer_stop();
 
       std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
+      auto name = as[level].actions.at(a)->action_name();
       if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
+
       force = FieldImplementation::projectForce(force); // Ta for gauge fields
       double end_force = usecond();
-      Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites());
+
-      std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] Force average: " << force_abs << std::endl;
+      MomFilter->applyFilter(force);
+      std::cout << GridLogIntegrator << " update_P : Level [" << level <<"]["<<a <<"] "<<name<< std::endl;
+      DumpSliceNorm("force ",force,Nd-1);
+      
+      Real force_abs   = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm.  nb. norm2(latt) = \sum_x norm2(latt[x]) 
+      Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR;    
+
+      Real force_max   = std::sqrt(maxLocalNorm2(force));
+      Real impulse_max = force_max * ep * HMC_MOMENTUM_DENOMINATOR;    
+
+      as[level].actions.at(a)->deriv_log(force_abs,force_max);
+      
+      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force average: " << force_abs <<" "<<name<<std::endl;
+      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force max    : " << force_max <<" "<<name<<std::endl;
+      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt average  : " << impulse_abs <<" "<<name<<std::endl;
+      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt max      : " << impulse_max <<" "<<name<<std::endl;
+
       Mom -= force * ep* HMC_MOMENTUM_DENOMINATOR;; 
       double end_full = usecond();
       double time_full  = (end_full - start_full) / 1e3;
       double time_force = (end_force - start_force) / 1e3;
       std::cout << GridLogMessage << "["<<level<<"]["<<a<<"] P update elapsed time: " << time_full << " ms (force: " << time_force << " ms)"  << std::endl;
+
     }
 
     // Force from the other representations
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
 
-    MomFilter->applyFilter(Mom);
   }
 
   void update_U(Field& U, double ep) 
@@ -161,8 +184,12 @@ protected:
   
   void update_U(MomentaField& Mom, Field& U, double ep) 
   {
+    MomentaField MomFiltered(Mom.Grid());
+    MomFiltered = Mom;
+    MomFilter->applyFilter(MomFiltered);
+
     // exponential of Mom*U in the gauge fields case
-    FieldImplementation::update_field(Mom, U, ep);
+    FieldImplementation::update_field(MomFiltered, U, ep);
 
     // Update the smeared fields, can be implemented as observer
     Smearer.set_Field(U);
@@ -205,6 +232,66 @@ public:
   const MomentaField & getMomentum() const{ return P; }
   
 
+  void reset_timer(void)
+  {
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+        as[level].actions.at(actionID)->reset_timer();
+      }
+    }
+  }
+  void print_timer(void)
+  {
+    std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
+    std::cout << GridLogMessage << " Refresh cumulative timings "<<std::endl;
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+	std::cout << GridLogMessage 
+		  << as[level].actions.at(actionID)->action_name()
+		  <<"["<<level<<"]["<< actionID<<"] "
+		  << as[level].actions.at(actionID)->refresh_us*1.0e-6<<" s"<< std::endl;
+      }
+    }
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    std::cout << GridLogMessage << " Action cumulative timings "<<std::endl;
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+	std::cout << GridLogMessage 
+		  << as[level].actions.at(actionID)->action_name()
+		  <<"["<<level<<"]["<< actionID<<"] "
+		  << as[level].actions.at(actionID)->S_us*1.0e-6<<" s"<< std::endl;
+      }
+    }
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    std::cout << GridLogMessage << " Force cumulative timings "<<std::endl;
+    std::cout << GridLogMessage << "------------------------- "<<std::endl;
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+	std::cout << GridLogMessage 
+		  << as[level].actions.at(actionID)->action_name()
+		  <<"["<<level<<"]["<< actionID<<"] "
+		  << as[level].actions.at(actionID)->deriv_us*1.0e-6<<" s"<< std::endl;
+      }
+    }
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    std::cout << GridLogMessage << " Force average size "<<std::endl;
+    std::cout << GridLogMessage << "------------------------- "<<std::endl;
+    for (int level = 0; level < as.size(); ++level) {
+      for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+	std::cout << GridLogMessage 
+		  << as[level].actions.at(actionID)->action_name()
+		  <<"["<<level<<"]["<< actionID<<"] : "
+		  <<" force max " << as[level].actions.at(actionID)->deriv_max_average()
+		  <<" norm "      << as[level].actions.at(actionID)->deriv_norm_average()
+		  <<" calls "     << as[level].actions.at(actionID)->deriv_num
+		  << std::endl;
+      }
+    }
+    std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
+  }
+  
   void print_parameters()
   {
     std::cout << GridLogMessage << "[Integrator] Name : "<< integrator_name() << std::endl;
@@ -223,7 +310,6 @@ public:
       }
     }
     std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
-
   }
 
   void reverse_momenta()
@@ -266,15 +352,19 @@ public:
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
         // get gauge field from the SmearingPolicy and
         // based on the boolean is_smeared in actionID
+	auto name = as[level].actions.at(actionID)->action_name();
+        std::cout << GridLogMessage << "refresh [" << level << "][" << actionID << "] "<<name << std::endl;
+
         Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+	as[level].actions.at(actionID)->refresh_timer_start();
         as[level].actions.at(actionID)->refresh(Us, sRNG, pRNG);
+	as[level].actions.at(actionID)->refresh_timer_stop();
       }
 
       // Refresh the higher representation actions
       as[level].apply(refresh_hireps, Representations, sRNG, pRNG);
     }
 
-    MomFilter->applyFilter(P);
   }
 
   // to be used by the actionlevel class to iterate
@@ -309,7 +399,9 @@ public:
         // based on the boolean is_smeared in actionID
         Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] action eval " << std::endl;
+	        as[level].actions.at(actionID)->S_timer_start();
         Hterm = as[level].actions.at(actionID)->S(Us);
+   	        as[level].actions.at(actionID)->S_timer_stop();
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] H = " << Hterm << std::endl;
         H += Hterm;
       }

Grid/stencil/Stencil.h

@@ -1120,7 +1120,7 @@ public:
 		    bytes);
 	}
 	
-	if ( compress.DecompressionStep() && comms_recv ) {
+	if ( compress.DecompressionStep()  ) {
 	  AddDecompress(&this->u_recv_buf_p[u_comm_offset],
 			&recv_buf[u_comm_offset],
 			words,Decompressions);
@@ -1206,8 +1206,8 @@ public:
 				  face_table[face_idx].size()*sizeof(face_table_host[0]));
 	}
 
-	if ( comms_send )
+	//	if ( comms_send )
-	  Gather_plane_exchange_table(face_table[face_idx],rhs,spointers,dimension,sx,cbmask,compress,permute_type);
+	Gather_plane_exchange_table(face_table[face_idx],rhs,spointers,dimension,sx,cbmask,compress,permute_type);
 	face_idx++;
 
 	//spointers[0] -- low
@@ -1239,7 +1239,7 @@ public:
 	    rpointers[i] = rp;
 
 	    int duplicate = CheckForDuplicate(dimension,sx,nbr_proc,(void *)rp,i,bytes,cbmask);
-	    if ( (!duplicate) ) { // Force comms for now
+	    if ( !duplicate  ) { 
 	      AddPacket((void *)sp,(void *)rp,
 			xmit_to_rank,comms_send,
 			recv_from_rank,comms_recv,
@@ -1253,9 +1253,7 @@ public:
 	  }
 	}
 
-	if ( comms_recv ) {
+	AddMerge(&this->u_recv_buf_p[u_comm_offset],rpointers,reduced_buffer_size,permute_type,Mergers);
-	  AddMerge(&this->u_recv_buf_p[u_comm_offset],rpointers,reduced_buffer_size,permute_type,Mergers);
-	}
 
 	u_comm_offset     +=buffer_size;
       }

Grid/tensors/Tensor_exp.h

@@ -55,7 +55,7 @@ template<class vtype, int N> accelerator_inline iVector<vtype, N> Exponentiate(c
 
 
 // Specialisation: Cayley-Hamilton exponential for SU(3)
-#ifndef GRID_CUDA
+#ifndef GRID_ACCELERATED
 template<class vtype, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0>::type * =nullptr> 
 accelerator_inline iMatrix<vtype,3> Exponentiate(const iMatrix<vtype,3> &arg, RealD alpha  , Integer Nexp = DEFAULT_MAT_EXP )
 {

Grid/threads/Accelerator.h

@@ -441,7 +441,7 @@ inline void acceleratorMemSet(void *base,int value,size_t bytes) { hipMemset(bas
 
 inline void acceleratorCopyDeviceToDeviceAsynch(void *from,void *to,size_t bytes) // Asynch
 {
-  hipMemcpyAsync(to,from,bytes, hipMemcpyDeviceToDevice,copyStream);
+  hipMemcpy(to,from,bytes, hipMemcpyDeviceToDevice);
 }
 inline void acceleratorCopySynchronise(void) { hipStreamSynchronize(copyStream); };
 
@@ -461,6 +461,8 @@ inline void acceleratorCopySynchronise(void) { hipStreamSynchronize(copyStream);
   accelerator_for2dNB(iter1, num1, iter2, num2, nsimd, { __VA_ARGS__ } ); \
   accelerator_barrier(dummy);
 
+#define GRID_ACCELERATED
+
 #endif
 
 //////////////////////////////////////////////

benchmarks/Benchmark_dwf_fp32.cc

@@ -143,7 +143,7 @@ void Benchmark(int Ls, Coordinate Dirichlet)
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedUniqueString(std::string("The 5D RNG"));
 
   LatticeFermionF src   (FGrid); random(RNG5,src);
-#if 1
+#if 0
   src = Zero();
   {
     Coordinate origin({0,0,0,latt4[2]-1,0});

systems/Spock/config-command

@@ -6,6 +6,8 @@
 --enable-simd=GPU \
 --disable-fermion-reps \
 --disable-gparity \
+--with-gmp=$OLCF_GMP_ROOT \
+--with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
 CXX=hipcc MPICXX=mpicxx \
 CXXFLAGS="-fPIC -I/opt/rocm-4.3.0/include/ -std=c++14 -I${MPICH_DIR}/include " \
 --prefix=/ccs/home/chulwoo/Grid \

systems/Spock/dwf8.slurm

@@ -1,8 +1,7 @@
 #!/bin/bash
 # Begin LSF Directives
 #SBATCH -A LGT104
-#SBATCH -t 01:00:00
+#SBATCH -t 3:00:00
-##SBATCH -U openmpThu
 #SBATCH -p ecp
 #SBATCH -J DWF
 #SBATCH -o DWF.%J
@@ -14,13 +13,12 @@ DIR=.
 module list
 export MPIR_CVAR_GPU_EAGER_DEVICE_MEM=0
 export MPICH_GPU_SUPPORT_ENABLED=1
-#export MPICH_SMP_SINGLE_COPY_MODE=XPMEM
+export MPICH_SMP_SINGLE_COPY_MODE=CMA
-export MPICH_SMP_SINGLE_COPY_MODE=NONE
+
-#export MPICH_SMP_SINGLE_COPY_MODE=CMA
 export OMP_NUM_THREADS=8
 
 AT=8
 echo MPICH_SMP_SINGLE_COPY_MODE $MPICH_SMP_SINGLE_COPY_MODE
-PARAMS=" --accelerator-threads ${AT} --grid 32.64.64.64 --mpi 1.2.2.2 --comms-overlap --shm 2048 --shm-mpi 0"
+PARAMS=" --accelerator-threads ${AT} --grid 16.16.16.48 --mpi 1.2.2.2 --comms-overlap --shm 2048 --shm-mpi 0"
-srun -n8 --label -c$OMP_NUM_THREADS --gpus-per-task=1 ./mpiwrapper.sh ./benchmarks/Benchmark_dwf_fp32 $PARAMS
+srun -N2 -n8 --label -c$OMP_NUM_THREADS --gpus-per-task=1 ./mpiwrapper.sh ./HMC/Mobius2p1f_DD_RHMC $PARAMS
 

systems/Spock/sourceme.sh

@@ -1,5 +1,9 @@
+module load emacs
 module load PrgEnv-gnu
 module load rocm/4.5.0
 module load gmp
 module load cray-fftw
 module load craype-accel-amd-gfx908
+export MPIR_CVAR_GPU_EAGER_DEVICE_MEM=0
+export MPICH_GPU_SUPPORT_ENABLED=1
+export LD_LIBRARY_PATH=/opt/cray/pe/gcc/mpfr/3.1.4/lib/:$LD_LIBRARY_PATH