finer timers in Benchmark_IO

Grid/parallelIO/BinaryIO.h

@@ -271,7 +271,7 @@ class BinaryIO {
 			      uint32_t &scidac_csumb)
   {
     grid->Barrier();
-    GridStopWatch timer; 
+    GridStopWatch timer, insideTimer; 
     GridStopWatch bstimer;
     
     nersc_csum=0;
@@ -363,7 +363,10 @@ class BinaryIO {
 	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< std::endl;
 	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);    assert(ierr==0);
 	ierr=MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);    assert(ierr==0);
-	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);    assert(ierr==0);
+  insideTimer.Start();
+	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);
+  insideTimer.Stop();
+  assert(ierr==0);
 	MPI_File_close(&fh);
 	MPI_Type_free(&fileArray);
 	MPI_Type_free(&localArray);
@@ -438,7 +441,9 @@ class BinaryIO {
         assert(ierr == 0);
 
         std::cout << GridLogDebug << "MPI write I/O write all " << file << std::endl;
+        insideTimer.Start();
         ierr = MPI_File_write_all(fh, &iodata[0], 1, localArray, &status);
+        insideTimer.Stop();
         assert(ierr == 0);
 
         MPI_Offset os;
@@ -516,8 +521,13 @@ class BinaryIO {
     if ( control & BINARYIO_READ) std::cout << " read  ";
     else                          std::cout << " write ";
     uint64_t bytes = sizeof(fobj)*iodata.size()*nrank;
-    std::cout<< lastPerf.size <<" bytes in "<< timer.Elapsed() <<" "
+    std::cout<< lastPerf.size <<"bytes in "<< timer.Elapsed() <<" "
 	     << lastPerf.mbytesPerSecond <<" MB/s "<<std::endl;
+    std::cout << GridLogMessage << "IOobject: pure MPI IO call " 
+              << lastPerf.size <<" bytes in " 
+              << insideTimer.Elapsed() << " "
+              << lastPerf.size/1024./1024./(insideTimer.useconds()/1.0e6)
+              <<" MB/s "<<std::endl;
 
     std::cout<<GridLogMessage<<"IOobject: endian and checksum overhead "<<bstimer.Elapsed()  <<std::endl;
 

Grid/qcd/smearing/WilsonFlow.h

@@ -6,8 +6,7 @@ Source file: ./lib/qcd/modules/plaquette.h
 
 Copyright (C) 2017
 
-Author: Guido Cossu  <guido.cossu@ed.ac.uk>
-Author: Chulwoo Jung <chulwoo@bnl.gov>
+Author: Guido Cossu <guido.cossu@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
@@ -36,7 +35,7 @@ template <class Gimpl>
 class WilsonFlow: public Smear<Gimpl>{
   unsigned int Nstep;
   unsigned int measure_interval;
-  mutable RealD epsilon, taus,tolerance;
+  mutable RealD epsilon, taus;
 
 
   mutable WilsonGaugeAction<Gimpl> SG;
@@ -48,15 +47,13 @@ class WilsonFlow: public Smear<Gimpl>{
 public:
   INHERIT_GIMPL_TYPES(Gimpl)
 
-  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1, RealD tol = 1e-3):
+  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
   Nstep(Nstep),
     epsilon(epsilon),
-    tolerance(tol),
     measure_interval(interval),
     SG(WilsonGaugeAction<Gimpl>(3.0)) {
     // WilsonGaugeAction with beta 3.0
     assert(epsilon > 0.0);
-    assert(tolerance > 0.0);
     LogMessage();
   }
 
@@ -67,8 +64,6 @@ public:
 	      << "[WilsonFlow] epsilon : " << epsilon << std::endl;
     std::cout << GridLogMessage
 	      << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
-    std::cout << GridLogMessage
-	      << "[WilsonFlow] tolerance : " << tolerance << std::endl;
   }
 
   virtual void smear(GaugeField&, const GaugeField&) const;
@@ -111,14 +106,11 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   if (maxTau - taus < epsilon){
     epsilon = maxTau-taus;
   }
-  std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
+  //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());
-
+  GaugeField tmp(U.Grid()), Uprime(U.Grid());
   Uprime = U;
-  Usave = U;
-
   SG.deriv(U, Z);
   Zprime = -Z;
   Z *= 0.25;                                  // Z0 = 1/4 * F(U)
@@ -136,33 +128,18 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   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 arXiv:1301.4388
+  // 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;
-// Wrong
-//  RealD diff = norm2(diffU);
-//  std::cout << GridLogMessage << "norm2: " << diff << std::endl;
-  
-//  RealD tol=1e-3;
-  
-  RealD diff = real(rankInnerMax(diffU,diffU));
-  diff = sqrt(diff)/18.; // distance defined in Ramos 
-
-  GridBase *grid = diffU.Grid();
-  std::cout << GridLogMessage << "max: " << diff << std::endl;
-  grid->GlobalMax(diff);
-  std::cout << GridLogMessage << "max: " << diff << std::endl;
+  RealD diff = norm2(diffU);
+  // adjust integration step
     
-  if(diff < tolerance) {
   taus += epsilon;
-//  std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
-  } else {
-    U = Usave;
-  }
+  //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
     
-  epsilon = epsilon*0.95*std::pow(tolerance/diff,1./3.);
-  std::cout << GridLogMessage << "Distance : "<<diff<<"New epsilon : " << epsilon << std::endl;
+  epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
+  //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
 
 }
 
@@ -207,11 +184,8 @@ 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;
-  taus = 0.;
+  taus = epsilon;
   unsigned int step = 0;
-  double measTau = epsilon*measure_interval;
-  std::cout << GridLogMessage << "measTau :"<< measTau << std::endl;
   do{
     step++;
     //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
@@ -219,12 +193,10 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
     std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
 		  << step << "  " << taus << "  "
 	      << energyDensityPlaquette(out) << std::endl;
-//    if( step % measure_interval == 0){
-    if( taus >  measTau ) {
+    if( step % measure_interval == 0){
       std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
 		<< step << "  " 
 		<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
-      measTau += epsilon*measure_interval;
     }
   } while (taus < maxTau);
 

tests/smearing/Test_WilsonFlow.cc

@@ -33,7 +33,6 @@ namespace Grid{
     GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
             int, steps,
             double, step_size,
-            double, tol,
             int, meas_interval,
             double, maxTau); // for the adaptive algorithm
        
@@ -83,27 +82,13 @@ int main(int argc, char **argv) {
   SU<Nc>::HotConfiguration(pRNG, Umu);
   
   typedef Grid::XmlReader       Serialiser;
-//  Serialiser Reader("input.xml");
-//  WFParameters WFPar(Reader);
-//  ConfParameters CPar(Reader);
-//  WFParameters WFPar;
-  int steps = 800;
-  double step_size=0.02;
-  double tol=1e-4;
-  int meas_interval=50;
-  double maxTau = 16;
-//  ConfParameters CPar;
-//  CPar. conf_prefix="configurations/ckpoint_lat";
-//  CPar. rng_prefix="rngs/ckpoint_rng";
-//  CPar. StartConfiguration=100,
-//  CPar. EndConfiguration=110,
-//  CPar. Skip=1;
-//  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
-  CheckpointerParameters CPPar("configurations/ckpoint_lat","rngs/ckpoint_rng");
+  Serialiser Reader("input.xml");
+  WFParameters WFPar(Reader);
+  ConfParameters CPar(Reader);
+  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
   BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
 
-//  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
-  for (int conf = 100; conf <= 110; conf+= 1){
+  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
 
   CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
 
@@ -111,10 +96,9 @@ int main(int argc, char **argv) {
   std::cout << GridLogMessage << "Initial plaquette: "
     << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
 
-  WilsonFlow<PeriodicGimplR> WF(steps, step_size, meas_interval);
+  WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
 
-//  WF.smear_adaptive(Uflow, Umu, maxTau);
-  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);