From 339be37dbad38de7423ff2488719f01764fd0196 Mon Sep 17 00:00:00 2001 From: neo Date: Wed, 13 Apr 2016 17:00:14 +0900 Subject: [PATCH] Debugging smeared HMC --- .../action/pseudofermion/TwoFlavourEvenOdd.h | 2 +- lib/qcd/hmc/HmcRunner.h | 4 +- lib/qcd/hmc/integrators/Integrator.h | 4 +- lib/qcd/smearing/APEsmearing.h | 12 ++++-- lib/qcd/smearing/GaugeConfiguration.h | 39 ++++++++++++++----- lib/qcd/smearing/StoutSmearing.h | 13 +++++-- lib/qcd/utils/WilsonLoops.h | 22 +++++------ tests/Test_hmc_EOWilsonFermionGauge.cc | 2 +- 8 files changed, 63 insertions(+), 35 deletions(-) diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 601bf9ec..16b2d814 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -100,7 +100,7 @@ namespace Grid{ PhiOdd =PhiOdd*scale; PhiEven=PhiEven*scale; - + }; ////////////////////////////////////////////////////// diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h index cb3428b4..b30da24b 100644 --- a/lib/qcd/hmc/HmcRunner.h +++ b/lib/qcd/hmc/HmcRunner.h @@ -94,7 +94,7 @@ public: // Smearing policy std::cout << GridLogMessage << " Creating the Stout class\n"; double rho = 0.1; // smearing parameter - int Nsmear = 3; // number of smearing levels + int Nsmear = 1; // number of smearing levels Smear_Stout Stout(rho); std::cout << GridLogMessage << " Creating the SmearedConfiguration class\n"; SmearedConfiguration SmearingPolicy(UGrid, Nsmear, Stout); @@ -144,7 +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"; + 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 eb1eeb3a..e8950f22 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -115,10 +115,12 @@ namespace Grid{ void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){ for(int a=0; aderiv(U,force); // deriv should not include Ta + GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); + as[level].actions.at(a)->deriv(Us,force); // deriv should not include Ta std::cout<is_smeared <is_smeared) Smearer.smeared_force(force); force = Ta(force); + std::cout<gSites()) <::avgPlaquette(previous_u); + std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl; + } } @@ -73,14 +77,14 @@ namespace Grid { pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu); pokeLorentz(iLambda, iLambda_mu, mu); } - StoutSmearing.derivative(SigmaK, iLambda, GaugeK); + StoutSmearing.derivative(SigmaK, iLambda, GaugeK);// derivative of SmearBase return SigmaK; } /*! @brief Returns smeared configuration at level 'Level' */ const GaugeField& get_smeared_conf(int Level) const{ return SmearedSet[Level]; } - + //==================================================================== void set_iLambda(GaugeLinkField& iLambda, GaugeLinkField& e_iQ, const GaugeLinkField& iQ, @@ -117,8 +121,8 @@ namespace Grid { w2 = w * w; cosw = cos(w); - emiu = toComplex(cos(u)) - timesI(toComplex(u)); - e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u)); + emiu = toComplex(cos(u)) - timesI(toComplex(sin(u))); + e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(sin(2.0*u))); r01 = (toComplex(2.0*u) + timesI(toComplex(2.0*(u2-w2)))) * e2iu + emiu * (toComplex(16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) + @@ -174,7 +178,7 @@ namespace Grid { } - + //==================================================================== public: GaugeField* ThinLinks; /*!< @brief Pointer to the thin links configuration */ @@ -201,12 +205,14 @@ namespace Grid { // attach the smeared routines to the thin links U and fill the smeared set void set_GaugeField(GaugeField& U){ fill_smearedSet(U);} + //==================================================================== void smeared_force(GaugeField& SigmaTilde) const{ if (smearingLevels > 0){ - GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid + GaugeField force = SigmaTilde;//actually = U*SigmaTilde GaugeLinkField tmp_mu(SigmaTilde._grid); for (int mu = 0; mu < Nd; mu++){ + // to get SigmaTilde tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu); pokeLorentz(force, tmp_mu, mu); } @@ -221,7 +227,8 @@ namespace Grid { } }// if smearingLevels = 0 do nothing } - + //==================================================================== + GaugeField& get_SmearedU(){ return SmearedSet[smearingLevels-1]; @@ -230,10 +237,22 @@ namespace Grid { GaugeField& get_U(bool smeared=false) { // get the config, thin links by default if (smeared){ - if (smearingLevels) return get_SmearedU(); - else return *ThinLinks; + if (smearingLevels){ + RealD impl_plaq = WilsonLoops::avgPlaquette(SmearedSet[smearingLevels-1]); + std::cout<< GridLogDebug << "getting U Plaq: " << impl_plaq<< std::endl; + return get_SmearedU(); + + } + else { + RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); + std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; + return *ThinLinks; + } } - else return *ThinLinks; + else{ + RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); + std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; + return *ThinLinks;} } }; diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index ede204c3..5b5e4ba3 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -33,7 +33,7 @@ namespace Grid { void smear(GaugeField& u_smr,const GaugeField& U) const{ GaugeField C(U._grid); - GaugeLinkField tmp(U._grid), q_mu(U._grid), Umu(U._grid); + GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid); std::cout<< GridLogDebug << "Stout smearing started\n"; @@ -42,8 +42,8 @@ namespace Grid { for (int mu = 0; mu U(4,grid); - for(int d=0;d(Umu,d); - } staple = zero; - - if(nu != mu) { - + GridBase *grid = Umu._grid; + + std::vector U(4,grid); + for(int d=0;d(Umu,d); + } + // mu // ^ // |__> nu - + // __ // | // __| // - + staple+=Gimpl::ShiftStaple( Gimpl::CovShiftForward (U[nu],nu, Gimpl::CovShiftBackward(U[mu],mu, diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index 13828c51..d9634666 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -66,7 +66,7 @@ public: TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp,CG,CG); - Nf2.is_smeared=false; + Nf2.is_smeared=true; //Collect actions ActionLevel Level1(1);