Adding tests for WilsonFlow classes

bootstrap.sh

@@ -1,4 +1,4 @@
-#!/usr/bin/env bash
+]#!/usr/bin/env bash
 
 EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.2.9.tar.bz2'
 

configure.ac

@@ -419,6 +419,7 @@ AC_CONFIG_FILES(tests/hadrons/Makefile)
 AC_CONFIG_FILES(tests/hmc/Makefile)
 AC_CONFIG_FILES(tests/solver/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
+AC_CONFIG_FILES(tests/smearing/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
 AC_CONFIG_FILES(extras/Makefile)
 AC_CONFIG_FILES(extras/Hadrons/Makefile)

lib/Init.cc

@@ -1,6 +1,6 @@
     /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+    Grid physics library, www.github.com/paboyle/Grid
 
     Source file: ./lib/Init.cc
 
@@ -36,7 +36,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <stdint.h>
 #include <unistd.h>
 #include <sys/mman.h>
-#include <sys/stat.h> 
+#include <sys/stat.h>
 #include <sys/time.h>
 #include <signal.h>
 #include <iostream>
@@ -92,14 +92,14 @@ const std::vector<int> GridDefaultSimd(int dims,int nsimd)
       if ( nn>=2) {
 	layout[d]=2;
 	nn/=2;
-      } else { 
+      } else {
 	layout[d]=1;
       }
     }
     assert(nn==1);
     return layout;
 }
-  
+
 ////////////////////////////////////////////////////////////
 // Command line parsing assist for stock controls
 ////////////////////////////////////////////////////////////
@@ -143,7 +143,7 @@ void GridCmdOptionIntVector(std::string &str,std::vector<int> & vec)
     vec.push_back(i);
     if(std::ispunct(ss.peek()))
       ss.ignore();
-  }    
+  }
   return;
 }
 
@@ -248,24 +248,24 @@ void Grid_init(int *argc,char ***argv)
     std::cout<<GridLogMessage<<"  --help : this message"<<std::endl;
     std::cout<<GridLogMessage<<std::endl;
     std::cout<<GridLogMessage<<"Geometry:"<<std::endl;
-    std::cout<<GridLogMessage<<"  --mpi n.n.n.n   : default MPI decomposition"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --mpi n.n.n.n   : default MPI decomposition"<<std::endl;
     std::cout<<GridLogMessage<<"  --threads n     : default number of OMP threads"<<std::endl;
-    std::cout<<GridLogMessage<<"  --grid n.n.n.n  : default Grid size"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --grid n.n.n.n  : default Grid size"<<std::endl;
-    std::cout<<GridLogMessage<<"  --shm  M        : allocate M megabytes of shared memory for comms"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --shm  M        : allocate M megabytes of shared memory for comms"<<std::endl;
     std::cout<<GridLogMessage<<std::endl;
     std::cout<<GridLogMessage<<"Verbose and debug:"<<std::endl;
     std::cout<<GridLogMessage<<"  --log list      : comma separted list of streams from Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours"<<std::endl;
-    std::cout<<GridLogMessage<<"  --decomposition : report on default omp,mpi and simd decomposition"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --decomposition : report on default omp,mpi and simd decomposition"<<std::endl;
     std::cout<<GridLogMessage<<"  --debug-signals : catch sigsegv and print a blame report"<<std::endl;
-    std::cout<<GridLogMessage<<"  --debug-stdout  : print stdout from EVERY node"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --debug-stdout  : print stdout from EVERY node"<<std::endl;
-    std::cout<<GridLogMessage<<"  --notimestamp   : suppress millisecond resolution stamps"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --notimestamp   : suppress millisecond resolution stamps"<<std::endl;
     std::cout<<GridLogMessage<<std::endl;
     std::cout<<GridLogMessage<<"Performance:"<<std::endl;
-    std::cout<<GridLogMessage<<"  --dslash-generic: Wilson kernel for generic Nc"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --dslash-generic: Wilson kernel for generic Nc"<<std::endl;
-    std::cout<<GridLogMessage<<"  --dslash-unroll : Wilson kernel for Nc=3"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --dslash-unroll : Wilson kernel for Nc=3"<<std::endl;
-    std::cout<<GridLogMessage<<"  --dslash-asm    : Wilson kernel for AVX512"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --dslash-asm    : Wilson kernel for AVX512"<<std::endl;
-    std::cout<<GridLogMessage<<"  --lebesgue      : Cache oblivious Lebesgue curve/Morton order/Z-graph stencil looping"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --lebesgue      : Cache oblivious Lebesgue curve/Morton order/Z-graph stencil looping"<<std::endl;
-    std::cout<<GridLogMessage<<"  --cacheblocking n.m.o.p : Hypercuboidal cache blocking"<<std::endl;    
+    std::cout<<GridLogMessage<<"  --cacheblocking n.m.o.p : Hypercuboidal cache blocking"<<std::endl;
     std::cout<<GridLogMessage<<std::endl;
     exit(EXIT_SUCCESS);
   }
@@ -281,22 +281,22 @@ void Grid_init(int *argc,char ***argv)
   std::string COL_BLUE   = GridLogColours.colour["BLUE"];
   std::string COL_YELLOW = GridLogColours.colour["YELLOW"];
   std::string COL_BACKGROUND = GridLogColours.colour["NORMAL"];
-  
+
   std::cout <<std::endl;
-  std::cout <<COL_RED  << "__|__|__|__|__"<<             "|__|__|_"<<COL_PURPLE<<"_|__|__|"<<                "__|__|__|__|__"<<std::endl; 
+  std::cout <<COL_RED  << "__|__|__|__|__"<<             "|__|__|_"<<COL_PURPLE<<"_|__|__|"<<                "__|__|__|__|__"<<std::endl;
-  std::cout <<COL_RED  << "__|__|__|__|__"<<             "|__|__|_"<<COL_PURPLE<<"_|__|__|"<<                "__|__|__|__|__"<<std::endl; 
+  std::cout <<COL_RED  << "__|__|__|__|__"<<             "|__|__|_"<<COL_PURPLE<<"_|__|__|"<<                "__|__|__|__|__"<<std::endl;
-  std::cout <<COL_RED  << "__|_ |  |  |  "<<             "|  |  | "<<COL_PURPLE<<" |  |  |"<<                "  |  |  | _|__"<<std::endl; 
+  std::cout <<COL_RED  << "__|_ |  |  |  "<<             "|  |  | "<<COL_PURPLE<<" |  |  |"<<                "  |  |  | _|__"<<std::endl;
-  std::cout <<COL_RED  << "__|_          "<<             "        "<<COL_PURPLE<<"        "<<                "          _|__"<<std::endl; 
+  std::cout <<COL_RED  << "__|_          "<<             "        "<<COL_PURPLE<<"        "<<                "          _|__"<<std::endl;
   std::cout <<COL_RED  << "__|_  "<<COL_GREEN<<" GGGG   "<<COL_RED<<" RRRR   "<<COL_BLUE  <<" III    "<<COL_PURPLE<<"DDDD  "<<COL_PURPLE<<"    _|__"<<std::endl;
   std::cout <<COL_RED  << "__|_  "<<COL_GREEN<<"G       "<<COL_RED<<" R   R  "<<COL_BLUE  <<"  I     "<<COL_PURPLE<<"D   D "<<COL_PURPLE<<"    _|__"<<std::endl;
   std::cout <<COL_RED  << "__|_  "<<COL_GREEN<<"G       "<<COL_RED<<" R   R  "<<COL_BLUE  <<"  I     "<<COL_PURPLE<<"D    D"<<COL_PURPLE<<"    _|__"<<std::endl;
   std::cout <<COL_BLUE << "__|_  "<<COL_GREEN<<"G  GG   "<<COL_RED<<" RRRR   "<<COL_BLUE  <<"  I     "<<COL_PURPLE<<"D    D"<<COL_GREEN <<"    _|__"<<std::endl;
   std::cout <<COL_BLUE << "__|_  "<<COL_GREEN<<"G   G   "<<COL_RED<<" R  R   "<<COL_BLUE  <<"  I     "<<COL_PURPLE<<"D   D "<<COL_GREEN <<"    _|__"<<std::endl;
   std::cout <<COL_BLUE << "__|_  "<<COL_GREEN<<" GGGG   "<<COL_RED<<" R   R  "<<COL_BLUE  <<" III    "<<COL_PURPLE<<"DDDD  "<<COL_GREEN <<"    _|__"<<std::endl;
-  std::cout <<COL_BLUE << "__|_          "<<             "        "<<COL_GREEN <<"        "<<                "          _|__"<<std::endl; 
+  std::cout <<COL_BLUE << "__|_          "<<             "        "<<COL_GREEN <<"        "<<                "          _|__"<<std::endl;
-  std::cout <<COL_BLUE << "__|__|__|__|__"<<             "|__|__|_"<<COL_GREEN <<"_|__|__|"<<                "__|__|__|__|__"<<std::endl; 
+  std::cout <<COL_BLUE << "__|__|__|__|__"<<             "|__|__|_"<<COL_GREEN <<"_|__|__|"<<                "__|__|__|__|__"<<std::endl;
-  std::cout <<COL_BLUE << "__|__|__|__|__"<<             "|__|__|_"<<COL_GREEN <<"_|__|__|"<<                "__|__|__|__|__"<<std::endl; 
+  std::cout <<COL_BLUE << "__|__|__|__|__"<<             "|__|__|_"<<COL_GREEN <<"_|__|__|"<<                "__|__|__|__|__"<<std::endl;
-  std::cout <<COL_BLUE << "  |  |  |  |  "<<             "|  |  | "<<COL_GREEN <<" |  |  |"<<                "  |  |  |  |  "<<std::endl; 
+  std::cout <<COL_BLUE << "  |  |  |  |  "<<             "|  |  | "<<COL_GREEN <<" |  |  |"<<                "  |  |  |  |  "<<std::endl;
   std::cout << std::endl;
   std::cout << std::endl;
   std::cout <<COL_YELLOW<< std::endl;
@@ -342,7 +342,7 @@ void Grid_init(int *argc,char ***argv)
   }
   if( GridCmdOptionExists(*argv,*argv+*argc,"--notimestamp") ){
     GridLogTimestamp(0);
-  } else { 
+  } else {
     GridLogTimestamp(1);
   }
 
@@ -350,7 +350,7 @@ void Grid_init(int *argc,char ***argv)
 		  Grid_default_latt,
 		  Grid_default_mpi);
 
-  std::cout << GridLogMessage << "Requesting "<< CartesianCommunicator::MAX_MPI_SHM_BYTES <<" byte stencil comms buffers "<<std::endl;
+  std::cout << GridLogDebug << "Requesting "<< CartesianCommunicator::MAX_MPI_SHM_BYTES <<" byte stencil comms buffers "<<std::endl;
 
   if( GridCmdOptionExists(*argv,*argv+*argc,"--decomposition") ){
     std::cout<<GridLogMessage<<"Grid Decomposition\n";
@@ -366,15 +366,22 @@ void Grid_init(int *argc,char ***argv)
   Grid_is_initialised = 1;
 }
 
-  
+
 void Grid_finalize(void)
 {
-#if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPI3) 
+#if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPI3)
   MPI_Finalize();
   Grid_unquiesce_nodes();
 #endif
 }
 
+void GridLogLayout() {
+    std::cout << GridLogMessage << "Grid Layout\n";
+    std::cout << GridLogMessage << "\tGlobal lattice size  : "<< GridCmdVectorIntToString(GridDefaultLatt()) << std::endl;
+    std::cout << GridLogMessage << "\tOpenMP threads       : "<< GridThread::GetThreads() <<std::endl;
+    std::cout << GridLogMessage << "\tMPI tasks            : "<< GridCmdVectorIntToString(GridDefaultMpi()) << std::endl;
+}
+
 void * Grid_backtrace_buffer[_NBACKTRACE];
 
 void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr)

lib/Init.h

@@ -1,6 +1,6 @@
     /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+    Grid physics library, www.github.com/paboyle/Grid
 
     Source file: ./lib/Init.h
 
@@ -46,6 +46,7 @@ namespace Grid {
   const int              &GridThreads(void)  ;
   void                    GridSetThreads(int t) ;
   void GridLogTimestamp(int);
+  void GridLogLayout();
 
   // Common parsing chores
   std::string GridCmdOptionPayload(char ** begin, char ** end, const std::string & option);

lib/qcd/smearing/WilsonFlow.h

@@ -38,9 +38,9 @@ class WilsonFlow: public Smear<Gimpl>{
     unsigned int Nstep;
     RealD epsilon;
 
-    WilsonGaugeAction<Gimpl> &SG;
+    mutable WilsonGaugeAction<Gimpl> SG;
 
-    void evolve_step(typename Gimpl::GaugeField&);
+    void evolve_step(typename Gimpl::GaugeField&) const;
     RealD tau(unsigned int t)const {return epsilon*(t+1.0); }
 
 
@@ -50,7 +50,7 @@ class WilsonFlow: public Smear<Gimpl>{
     explicit WilsonFlow(unsigned int Nstep, RealD epsilon):
         Nstep(Nstep),
         epsilon(epsilon),
-        SG(3.0) {
+        SG(WilsonGaugeAction<Gimpl>(3.0)) {
             // WilsonGaugeAction with beta 3.0
             assert(epsilon > 0.0);
             LogMessage();
@@ -58,24 +58,29 @@ class WilsonFlow: public Smear<Gimpl>{
 
     void LogMessage() {
         std::cout << GridLogMessage
-            << "[WilsonFlow] Nstep   : " << Nstep
+            << "[WilsonFlow] Nstep   : " << Nstep << std::endl;
-            << "[WilsonFlow] epsilon : " << epsilon
+        std::cout << GridLogMessage
-            << "[WilsonFlow] full trajectory : " << Nstep * epsilon
+            << "[WilsonFlow] epsilon : " << epsilon << std::endl;
-            << 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 {
-                  assert(0);
+        assert(0);
-                  // undefined for WilsonFlow
+        // undefined for WilsonFlow
     }
 
     RealD energyDensityPlaquette(unsigned int step, const GaugeField& U) const;
 };
 
+
+////////////////////////////////////////////////////////////////////////////////
+// Implementations
+////////////////////////////////////////////////////////////////////////////////
 template <class Gimpl>
-void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U) {
+void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U) const{
     GaugeField Z(U._grid);
     GaugeField tmp(U._grid);
     SG.deriv(U, Z);
@@ -104,13 +109,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
     out = in;
     for (unsigned int step = 0; step < Nstep; step++) {
         evolve_step(out);
-        // Energy density, plaquette
+        std::cout << "[WilsonFlow] Energy density (plaq) : "
-        std::cout << GridLogMessage
+            << step << " "
-            << "[WilsonFlow] Energy Density (step) : "
+            << energyDensityPlaquette(step,out) << std::endl;
-            << step
+    }
-            << " " << energyDensityPlaquette(step, out)
+
-            << std::endl;
-}
 }
 
 }  // namespace QCD

lib/qcd/utils/WilsonLoops.h

@@ -99,6 +99,7 @@ public:
   //////////////////////////////////////////////////
   static RealD sumPlaquette(const GaugeLorentz &Umu) {
     std::vector<GaugeMat> U(Nd, Umu._grid);
+    // inefficient here
     for (int mu = 0; mu < Nd; mu++) {
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }

tests/Makefile.am

@@ -1,4 +1,4 @@
-SUBDIRS = core forces hmc solver debug	
+SUBDIRS = core forces hmc solver debug smearing
 
 if BUILD_CHROMA_REGRESSION
   SUBDIRS+= qdpxx

tests/hmc/Test_hmc_WilsonGauge.cc

@@ -35,9 +35,7 @@ int main(int argc, char **argv) {
   using namespace Grid::QCD;
 
   Grid_init(&argc, &argv);
-  int threads = GridThread::GetThreads();
+  GridLogLayout();
-  // here make a routine to print all the relevant information on the run
-  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
 
    // Typedefs to simplify notation
   typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm

tests/smearing/Makefile.am

@@ -0,0 +1 @@
+include Make.inc

tests/smearing/Test_WilsonFlow.cc

@@ -0,0 +1,78 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/hmc/Test_WilsonFlow.cc
+
+Copyright (C) 2017
+
+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
+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>
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size, simd_layout, mpi_layout);
+
+  std::vector<int> seeds({1, 2, 3, 4, 5});
+  GridParallelRNG pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+
+  LatticeGaugeField Umu(&Grid), Uflow(&Grid);
+  SU<Nc>::HotConfiguration(pRNG, Umu);
+
+  std::cout << std::setprecision(15);
+  std::cout << GridLogMessage << "Plaquette: "
+    << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
+
+  WilsonFlow<PeriodicGimplR> WF(200, 0.01);
+
+  WF.smear(Uflow, Umu);
+
+  RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
+  std::cout << GridLogMessage << "Plaquette: "<< WFlow_plaq << std::endl;
+
+  std::cout<< GridLogMessage << " Admissibility check:\n";
+  const double sp_adm = 0.067;                // admissible threshold
+  const double pl_adm = 1.0-sp_adm/Nc;
+  std::cout << GridLogMessage << "   (pl_adm =" << pl_adm << ")\n";
+
+  // Need min and reduce min for this function
+  //double sp_max = NC_*(1.0-stpl.plaq_min(U,pl_adm));
+  double sp_ave = Nc*(1.0-WFlow_plaq);
+
+  //std::cout<< GridLogMessage << "   sp_max = "        << sp_max <<"\n";
+  std::cout<< GridLogMessage << "   sp_ave = "        << sp_ave <<"\n";
+  std::cout<< GridLogMessage << "   (sp_admissible = "<< sp_adm <<")\n";
+  //std::cout<< GridLogMessage << "   sp_admissible - sp_max = "<<sp_adm-sp_max <<"\n";
+  std::cout<< GridLogMessage << "   sp_admissible - sp_ave = "<<sp_adm-sp_ave <<"\n";
+
+  Grid_finalize();
+}  // main