Improved logging output for RHMC bounds checks

In GenericHMCRunner, exposed functionality for initializing gauge fields and RNG for external use

Grid/qcd/action/pseudofermion/Bounds.h

@@ -40,10 +40,11 @@ NAMESPACE_BEGIN(Grid);
       X=X-Y;
       RealD Nd = norm2(X);
       std::cout << "************************* "<<std::endl;
-      std::cout << " noise                         = "<<Nx<<std::endl;
+      std::cout << " | noise |^2                         = "<<Nx<<std::endl;
-      std::cout << " (MdagM^-1/2)^2  noise         = "<<Nz<<std::endl;
+      std::cout << " | (MdagM^-1/2)^2  noise |^2         = "<<Nz<<std::endl;
-      std::cout << " MdagM (MdagM^-1/2)^2  noise   = "<<Ny<<std::endl;
+      std::cout << " | MdagM (MdagM^-1/2)^2  noise |^2   = "<<Ny<<std::endl;
-      std::cout << " noise - MdagM (MdagM^-1/2)^2  noise   = "<<Nd<<std::endl;
+      std::cout << " | noise - MdagM (MdagM^-1/2)^2  noise |^2  = "<<Nd<<std::endl;
+      std::cout << " | noise - MdagM (MdagM^-1/2)^2  noise|/|noise| = " << std::sqrt(Nd/Nx) << std::endl;
       std::cout << "************************* "<<std::endl;
       assert( (std::sqrt(Nd/Nx)<tol) && " InverseSqrtBoundsCheck ");
     }
@@ -90,10 +91,11 @@ NAMESPACE_BEGIN(Grid);
       X=X-Y;
       RealD Nd = norm2(X);
       std::cout << "************************* "<<std::endl;
-      std::cout << " noise                         = "<<Nx<<std::endl;
+      std::cout << " | noise |^2                         = "<<Nx<<std::endl;
-      std::cout << " (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise         = "<<Nz<<std::endl;
+      std::cout << " | (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2        = "<<Nz<<std::endl;
-      std::cout << " MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise   = "<<Ny<<std::endl;
+      std::cout << " | MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2   = "<<Ny<<std::endl;
-      std::cout << " noise - MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise   = "<<Nd<<std::endl;
+      std::cout << " | noise - MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2  = "<<Nd<<std::endl;
+      std::cout << " | noise - MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |/| noise |  = "<<std::sqrt(Nd/Nx)<<std::endl;
       std::cout << "************************* "<<std::endl;
       assert( (std::sqrt(Nd/Nx)<tol) && " InversePowerBoundsCheck ");
     }

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h

@@ -180,6 +180,8 @@ NAMESPACE_BEGIN(Grid);
 	return sstream.str();
       }
 
+      //Access the fermion field
+      const FermionField &getPhiOdd() const{ return PhiOdd; }
       
       virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
@@ -189,21 +191,19 @@ NAMESPACE_BEGIN(Grid);
 	//        = e^{- phi^dag  (VdagV)^1/(2*inv_pow) (MdagM)^-1/(2*inv_pow) (MdagM)^-1/(2*inv_pow)  (VdagV)^1/(2*inv_pow) phi}
 	//
 	// Phi =  (VdagV)^-1/(2*inv_pow) Mdag^{1/(2*inv_pow)} eta 
-	//
-	// P(eta) = e^{- eta^dag eta}
-	//	
-	// General gaussian random draws from   e^{x^2/(2 sig^2)} => sig^2 = 0.5.
-	// 
-	// So eta should be of width sig = 1/sqrt(2).
 	
 	std::cout<<GridLogMessage << action_name() << " refresh: starting" << std::endl;
-	RealD scale = std::sqrt(0.5);
 
 	FermionField eta(NumOp.FermionGrid());
 	FermionField etaOdd (NumOp.FermionRedBlackGrid());
 	FermionField etaEven(NumOp.FermionRedBlackGrid());
 	FermionField     tmp(NumOp.FermionRedBlackGrid());
 
+	// P(eta) \propto e^{- eta^dag eta}
+	//	
+	// The gaussian function draws from  P(x) \propto e^{- x^2 / 2 }    [i.e. sigma=1]
+	// Thus eta = x/sqrt{2} = x * sqrt(1/2)
+	RealD scale = std::sqrt(0.5);
 	gaussian(pRNG,eta);	eta=eta*scale;
 
 	pickCheckerboard(Even,etaEven,eta);

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,22 @@ 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);
 
     Smearing.set_Field(U);
 

Grid/qcd/hmc/HMCResourceManager.h

@@ -227,6 +227,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
     // the default is the first grid registered