Debugging Smearing routines (set_fj)

lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -189,7 +189,7 @@ namespace Grid{
 	assert(DenOp.ConstEE() == 1);
 
 	//dSdU = -Ta(dSdU);
-	dSDu = -dSdU;
+	dSdU = -dSdU;
 	
       };
     };

lib/qcd/hmc/HmcRunner.h

@@ -91,7 +91,15 @@ public:
 
     
     // Create integrator, including the smearing policy
-    SmearedConfiguration<Gimpl> SmearingPolicy; // simplest empty smearer, construct here more complex smearers
+    // Smearing policy
+    std::cout << GridLogMessage << " Creating the Stout class\n";
+    double rho = 0.1; // smearing parameter
+    int Nsmear = 3;   // number of smearing levels
+    Smear_Stout<Gimpl> Stout(rho);
+    std::cout << GridLogMessage << " Creating the SmearedConfiguration class\n";
+    SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
+    std::cout << GridLogMessage << " done\n";
+    //////////////
     typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl> >  IntegratorType;// change here to change the algorithm
     IntegratorParameters MDpar(20);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
@@ -136,6 +144,7 @@ public:
 
     // Attach the gauge field to the smearing Policy and create the fill the smeared set
     // notice that the unit configuration is singular in this procedure
+    std::cout << GridLogMessage << " Filling the smeared set\n";
     SmearingPolicy.set_GaugeField(U);
     
     HybridMonteCarlo<GaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U); 

lib/qcd/hmc/integrators/Integrator.h

@@ -116,6 +116,8 @@ namespace Grid{
 	for(int a=0; a<as[level].actions.size(); ++a){
 	  GaugeField force(U._grid);
 	  as[level].actions.at(a)->deriv(U,force); // deriv should not include Ta
+	  std::cout<<GridLogMessage<< "P Level: "<< level<< " action: "<<a<< " Smeared: "
+		   <<as[level].actions.at(a)->is_smeared <<std::endl;
 	  if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
 	  force = Ta(force);
 	  Mom = Mom - force*ep;
@@ -174,6 +176,9 @@ namespace Grid{
 	  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
+	    std::cout<<GridLogMessage<< "Refresh Level: "<< level<< " action: "<<actionID<< " Smeared: "
+		     <<as[level].actions.at(actionID)->is_smeared <<std::endl;
+
 	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
 	    as[level].actions.at(actionID)->refresh(Us, pRNG);
 	  }
@@ -201,8 +206,8 @@ namespace Grid{
 	    // get gauge field from the SmearingPolicy and
 	    // based on the boolean is_smeared in actionID
 	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+	    std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
 	    Hterm = as[level].actions.at(actionID)->S(Us);
-	    std::cout<<GridLogMessage << "Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
 	    H += Hterm;
 	  }
 	}

lib/qcd/smearing/GaugeConfiguration.h

@@ -35,17 +35,21 @@ namespace Grid {
       //====================================================================
       void fill_smearedSet(GaugeField& U){
 	ThinLinks = &U; //attach the smearing routine to the field U
-	//check that the pointer is not null
+
+	//check the pointer is not null
 	if (ThinLinks==NULL) 
 	  std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n";
 	
 	if (smearingLevels > 0){
-	  GaugeField previous_u(U._grid);
 	  std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
+	  GaugeField previous_u(ThinLinks->_grid);
+	  std::cout<< GridLogDebug << "[SmearedConfiguration] Loop\n";
 	  
 	  previous_u = *ThinLinks;
 	  for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
+	    std::cout<< GridLogDebug << "[SmearedConfiguration] Loop: "<< smearLvl << "\n";
 	    StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
+	    std::cout<< GridLogDebug << "[SmearedConfiguration] Loop assign: "<< smearLvl << "\n";
 	    previous_u = SmearedSet[smearLvl];
 	  }
 
@@ -201,22 +205,24 @@ namespace Grid {
       void set_GaugeField(GaugeField& U){ fill_smearedSet(U);}
       
       void smeared_force(GaugeField& SigmaTilde) const{
-	GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid
-	GaugeLinkField tmp_mu(SigmaTilde._grid);
-	
-	for (int mu = 0; mu < Nd; mu++){
-	  tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
+	if (smearingLevels > 0){
+	  GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid
+	  GaugeLinkField tmp_mu(SigmaTilde._grid);
+	  
+	  for (int mu = 0; mu < Nd; mu++){
+	    tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
 	    pokeLorentz(force, tmp_mu, mu);
-	}
-	for(int ismr = smearingLevels - 1; ismr > 0; --ismr)
-	  force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1));
-	
-	force = AnalyticSmearedForce(force,*ThinLinks);
-	
-	for (int mu = 0; mu < Nd; mu++){
-	  tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
-	  pokeLorentz(SigmaTilde, tmp_mu, mu);
-	}
+	  }
+	  for(int ismr = smearingLevels - 1; ismr > 0; --ismr)
+	    force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1));
+	  
+	  force = AnalyticSmearedForce(force,*ThinLinks);
+	  
+	  for (int mu = 0; mu < Nd; mu++){
+	    tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
+	    pokeLorentz(SigmaTilde, tmp_mu, mu);
+	  }
+	}// if smearingLevels = 0 do nothing
       }
       
       

lib/qcd/smearing/StoutSmearing.h

@@ -25,7 +25,7 @@ namespace Grid {
       }
       
       /*! Default constructor */
-      Smear_Stout():SmearBase(new Smear_APE < Gimpl > ()){
+      Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){
 	static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
       }
       
@@ -68,6 +68,7 @@ namespace Grid {
 
 	// only one Lorentz direction at a time 
 
+	std::cout<< GridLogDebug << "Stout smearing exponentiate iQ\n";
 	GridBase *grid = iQ._grid;
 	GaugeLinkField unity(grid);
 	unity=1.0;
@@ -92,7 +93,8 @@ namespace Grid {
 
 	GridBase *grid = u._grid;
 	LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
-	
+
+	std::cout<< GridLogDebug << "Stout smearing set uw\n";
 	
 	c0    = - toReal(imag(trace(iQ3))) * one_over_three;
 	c1    = - toReal(real(trace(iQ2))) * one_over_two;
@@ -107,6 +109,7 @@ namespace Grid {
       void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2,
 		  const LatticeReal& u, const LatticeReal& w) const{
 
+	std::cout<< GridLogDebug << "Stout smearing set fj\n";
 	GridBase *grid = u._grid;
 	LatticeReal xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid);
 	LatticeComplex fden(grid);
@@ -117,8 +120,10 @@ namespace Grid {
 	u2    = u * u;
 	w2    = w * w;
 	cosw  = cos(w);
-	
+
+	std::cout<< GridLogDebug << "Stout smearing first toComplex\n";
 	ixi0  = timesI(toComplex(xi0));
+	std::cout<< GridLogDebug << "Stout smearing others toComplex\n";
 	emiu  = toComplex(cos(u)) - timesI(toComplex(u));
 	e2iu  = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u));