diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h index b2ae4d21..edb6beaa 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h @@ -189,7 +189,7 @@ namespace Grid{ assert(DenOp.ConstEE() == 1); //dSdU = -Ta(dSdU); - dSDu = -dSdU; + dSdU = -dSdU; }; }; diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h index 6dac0d53..cb3428b4 100644 --- a/lib/qcd/hmc/HmcRunner.h +++ b/lib/qcd/hmc/HmcRunner.h @@ -91,7 +91,15 @@ public: // Create integrator, including the smearing policy - SmearedConfiguration 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 Stout(rho); + std::cout << GridLogMessage << " Creating the SmearedConfiguration class\n"; + SmearedConfiguration SmearingPolicy(UGrid, Nsmear, Stout); + std::cout << GridLogMessage << " done\n"; + ////////////// typedef MinimumNorm2 > 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 HMC(HMCpar, MDynamics,sRNG,pRNG,U); diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index 5a70a1ec..0a6fb9c7 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -116,6 +116,8 @@ namespace Grid{ for(int a=0; aderiv(U,force); // deriv should not include Ta + std::cout<is_smeared <is_smeared) Smearer.smeared_force(force); force = Ta(force); Mom = Mom - force*ep; @@ -174,6 +176,9 @@ namespace Grid{ for(int actionID=0; actionIDis_smeared <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< 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 } diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index dc4b29be..940dee1b 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/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)); diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index c7f2d534..d9634666 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -65,7 +65,9 @@ public: ConjugateGradient CG(1.0e-8,10000); TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp,CG,CG); - + + Nf2.is_smeared=true; + //Collect actions ActionLevel Level1(1); Level1.push_back(&Nf2);