Better logging and moving momentumfilter

Grid/qcd/hmc/HMC.h

@@ -116,22 +116,17 @@ private:
 
     random(sRNG, rn_test);
 
-    std::cout << GridLogMessage
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
-              << "--------------------------------------------------\n";
+    std::cout << GridLogMessage << "exp(-dH) = " << prob << "  Random = " << rn_test << "\n";
-    std::cout << GridLogMessage << "exp(-dH) = " << prob
+    std::cout << GridLogMessage << "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 << GridLogMessage
+      std::cout << GridLogMessage << "--------------------------------------------------\n";
-                << "--------------------------------------------------\n";
       return true;
     } else {  // rejected
       std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n";
-      std::cout << GridLogMessage
+      std::cout << GridLogMessage << "--------------------------------------------------\n";
-                << "--------------------------------------------------\n";
       return false;
     }
   }
@@ -140,18 +135,63 @@ 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";
+
+    //////////////////////////////////////////////////////////////////////////////////////////////////////
+    // 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.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){
@@ -164,18 +204,15 @@ private:
     }
     ///////////////////////////////////////////////////////////
 
-
     std::cout.precision(15);
-    std::cout << GridLogMessage << "Total H after trajectory  = " << H1
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
-	      << "  dH = " << H1 - H0 << "\n";
+    std::cout << GridLogMessage << "Total H after trajectory  = " << H1 << "  dH = " << H1 - H0 << "\n";
+    std::cout << GridLogMessage << "--------------------------------------------------\n";
     std::cout.precision(current_precision);
     
     return (H1 - H0);
   }
 
-
-  
-
 public:
   /////////////////////////////////////////
   // Constructor
@@ -196,8 +233,11 @@ 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";
+
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
       	std::cout << GridLogMessage << "-- Thermalization" << std::endl;
       }
@@ -217,8 +257,6 @@ public:
       if (accept)
         Ucur = Ucopy; 
       
-     
-      
       double t1=usecond();
       std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 

Grid/qcd/hmc/integrators/Integrator.h

@@ -33,7 +33,6 @@ directory
 #define INTEGRATOR_INCLUDED
 
 #include <memory>
-#include "MomentumFilter.h"
 
 NAMESPACE_BEGIN(Grid);
 

Grid/qcd/hmc/integrators/MomentumFilter.h

@@ -1,94 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h
-
-Copyright (C) 2015
-
-Author: Christopher Kelly <ckelly@bnl.gov>
-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 */
-//--------------------------------------------------------------------
-#ifndef MOMENTUM_FILTER
-#define MOMENTUM_FILTER
-
-NAMESPACE_BEGIN(Grid);
-
-//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh
-
-template<typename MomentaField>
-struct MomentumFilterBase{
-  virtual void applyFilter(MomentaField &P) const;
-};
-
-//Do nothing
-template<typename MomentaField>
-struct MomentumFilterNone: public MomentumFilterBase<MomentaField>{
-  void applyFilter(MomentaField &P) const override{}
-};
-
-//Multiply each site/direction by a Lorentz vector complex number field
-//Can be used to implement a mask, zeroing out sites
-template<typename MomentaField>
-struct MomentumFilterApplyPhase: public MomentumFilterBase<MomentaField>{
-  typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
-  typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
-  typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
-  typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
-  
-  LatticeLorentzScalarType phase;
- 
-  MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){}
-
-  //Default to uniform field of (1,0)
-  MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){
-    LorentzScalarType one;
-    for(int mu=0;mu<Nd;mu++)
-      one(mu)()() = scalar_type(1.);
-    
-    phase = one;
-  }
-
-  void applyFilter(MomentaField &P) const override{
-    conformable(P,phase);
-    autoView( P_v , P, AcceleratorWrite);
-    autoView( phase_v , phase, AcceleratorRead);
-
-    accelerator_for(ss,P_v.size(),MomentaField::vector_type::Nsimd(),{
-    	auto site_mom = P_v(ss);
-    	auto site_phase = phase_v(ss);
-	for(int mu=0;mu<Nd;mu++)
-	  site_mom(mu) = site_mom(mu) * site_phase(mu);
-    	coalescedWrite(P_v[ss], site_mom);
-      });
-    
-  }
-
-
-};
-
-
-
-
-NAMESPACE_END(Grid);
-
-#endif