From 092fa0d8dab94f5d177ac36ab6a36a07e80c1990 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 1 Jul 2016 16:06:20 +0100 Subject: [PATCH] Debugged set_fj, to be fixed: BUG in imag() --- lib/Log.h | 1 - lib/qcd/hmc/integrators/Integrator.h | 6 +- lib/qcd/smearing/APEsmearing.h | 116 +++++++-------- lib/qcd/smearing/GaugeConfiguration.h | 203 +++++++++++++------------- lib/qcd/smearing/StoutSmearing.h | 135 ++++++++--------- 5 files changed, 233 insertions(+), 228 deletions(-) diff --git a/lib/Log.h b/lib/Log.h index 9f721d33..c7fde889 100644 --- a/lib/Log.h +++ b/lib/Log.h @@ -63,7 +63,6 @@ namespace Grid { colour["WHITE"] ="\033[37m"; colour["NORMAL"] ="\033[0;39m"; } else { - std::cout << "Switching off colours\n"; colour["BLACK"] =""; colour["RED"] =""; colour["GREEN"] =""; diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index e8950f22..58e1f585 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -117,10 +117,10 @@ namespace Grid{ GaugeField force(U._grid); 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 <is_smeared) Smearer.smeared_force(force); force = Ta(force); - std::cout<gSites()) <gSites()) < - class Smear_APE: public Smear{ - private: + class Smear_APE: public Smear{ + private: const std::vector rho;/*!< Array of weights */ //This member must be private - we do not want to control from outside - std::vector set_rho(const double common_rho)const { - std::vector res; - - for(int mn=0; mn set_rho(const double common_rho)const { + std::vector res; - public: + for(int mn=0; mn& rho_):rho(rho_){} - Smear_APE(double rho_val):rho(set_rho(rho_val)){} - Smear_APE():rho(set_rho(1.0)){} - ~Smear_APE(){} + Smear_APE(const std::vector& rho_):rho(rho_){} + Smear_APE(double rho_val):rho(set_rho(rho_val)){} + Smear_APE():rho(set_rho(1.0)){} + ~Smear_APE(){} /////////////////////////////////////////////////////////////////////////////// - void smear(GaugeField& u_smr, const GaugeField& U)const{ - GridBase *grid = U._grid; - double d_rho; - GaugeLinkField Cup(grid), tmp_stpl(grid); - WilsonLoops WL; - u_smr = zero; + void smear(GaugeField& u_smr, const GaugeField& U)const{ + GridBase *grid = U._grid; + double d_rho; + GaugeLinkField Cup(grid), tmp_stpl(grid); + WilsonLoops WL; + u_smr = zero; - for(int mu=0; mugSites(); WilsonLoops WL; @@ -76,27 +78,27 @@ namespace Grid { GaugeLinkField U_mu(grid), U_nu(grid); GaugeLinkField sh_field(grid), temp_Sigma(grid); Real rho_munu, rho_numu; - - for(int mu = 0; mu < Nd; ++mu){ - U_mu = PeekIndex( U, mu); - iLambda_mu = PeekIndex(iLambda, mu); - - for(int nu = 0; nu < Nd; ++nu){ - if(nu==mu) continue; - U_nu = PeekIndex( U, nu); - iLambda_nu = PeekIndex(iLambda, nu); - - rho_munu = rho[mu + Nd * nu]; - rho_numu = rho[nu + Nd * mu]; - - WL.StapleUpper(staple, U, mu, nu); - temp_Sigma = -rho_numu*staple*iLambda_nu; + for(int mu = 0; mu < Nd; ++mu){ + U_mu = PeekIndex( U, mu); + iLambda_mu = PeekIndex(iLambda, mu); + + for(int nu = 0; nu < Nd; ++nu){ + if(nu==mu) continue; + U_nu = PeekIndex( U, nu); + iLambda_nu = PeekIndex(iLambda, nu); + + rho_munu = rho[mu + Nd * nu]; + rho_numu = rho[nu + Nd * mu]; + + WL.StapleUpper(staple, U, mu, nu); + + temp_Sigma = -rho_numu*staple*iLambda_nu; //-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x) - Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu); + Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu); sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity? - + temp_Sigma = rho_numu*sh_field*staple; //r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x) Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu); @@ -119,10 +121,10 @@ namespace Grid { sh_field = Cshift(temp_Sigma, nu, -1); Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu); - } } - } - }; +} +} +}; diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h index a87a7ef4..2202d545 100644 --- a/lib/qcd/smearing/GaugeConfiguration.h +++ b/lib/qcd/smearing/GaugeConfiguration.h @@ -6,10 +6,10 @@ #ifndef GAUGE_CONFIG_ #define GAUGE_CONFIG_ -namespace Grid { - - namespace QCD { - + namespace Grid { + + namespace QCD { + /*! @brief Smeared configuration container @@ -23,42 +23,43 @@ namespace Grid { */ template class SmearedConfiguration { - public: - INHERIT_GIMPL_TYPES(Gimpl) ; - - private: - const unsigned int smearingLevels; - Smear_Stout StoutSmearing; - std::vector SmearedSet; - + public: + INHERIT_GIMPL_TYPES(Gimpl) ; + + private: + const unsigned int smearingLevels; + Smear_Stout StoutSmearing; + std::vector SmearedSet; + // Member functions //==================================================================== - void fill_smearedSet(GaugeField& U){ + void fill_smearedSet(GaugeField& U){ ThinLinks = &U; //attach the smearing routine to the field U //check the pointer is not null if (ThinLinks==NULL) - std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n"; + std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n"; if (smearingLevels > 0){ - std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; - GaugeField previous_u(ThinLinks->_grid); - - previous_u = *ThinLinks; - for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){ - StoutSmearing.smear(SmearedSet[smearLvl],previous_u); - previous_u = SmearedSet[smearLvl]; + std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; + GaugeField previous_u(ThinLinks->_grid); - RealD impl_plaq = WilsonLoops::avgPlaquette(previous_u); - std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl; + previous_u = *ThinLinks; + for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){ + StoutSmearing.smear(SmearedSet[smearLvl],previous_u); + previous_u = SmearedSet[smearLvl]; - } + // For debug purposes + RealD impl_plaq = WilsonLoops::avgPlaquette(previous_u); + std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl; + + } } - } +} //==================================================================== - GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, - const GaugeField& GaugeK) const{ +GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, + const GaugeField& GaugeK) const{ GridBase *grid = GaugeK._grid; GaugeField C(grid), SigmaK(grid), iLambda(grid); GaugeLinkField iLambda_mu(grid); @@ -69,27 +70,27 @@ namespace Grid { StoutSmearing.BaseSmear(C, GaugeK); for (int mu = 0; mu < Nd; mu++){ - Cmu = peekLorentz( C,mu); - GaugeKmu = peekLorentz(GaugeK,mu); - SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu); - iQ = Ta(Cmu*adj(GaugeKmu)); - set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu); - pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu); - pokeLorentz(iLambda, iLambda_mu, mu); + Cmu = peekLorentz( C,mu); + GaugeKmu = peekLorentz(GaugeK,mu); + SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu); + iQ = Ta(Cmu*adj(GaugeKmu)); + set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu); + pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu); + pokeLorentz(iLambda, iLambda_mu, mu); } 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{ +const GaugeField& get_smeared_conf(int Level) const{ return SmearedSet[Level]; - } +} //==================================================================== - void set_iLambda(GaugeLinkField& iLambda, - GaugeLinkField& e_iQ, - const GaugeLinkField& iQ, - const GaugeLinkField& Sigmap, - const GaugeLinkField& GaugeK)const{ +void set_iLambda(GaugeLinkField& iLambda, + GaugeLinkField& e_iQ, + const GaugeLinkField& iQ, + const GaugeLinkField& Sigmap, + const GaugeLinkField& GaugeK)const{ GridBase *grid = iQ._grid; GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid); GaugeLinkField unity(grid); @@ -103,9 +104,9 @@ namespace Grid { LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid); LatticeComplex r22(grid), tr1(grid), tr2(grid); LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(grid); - LatticeReal unitReal(grid); + LatticeReal LatticeUnitReal(grid); - unitReal = 1.0; + LatticeUnitReal = 1.0; // Exponential iQ2 = iQ * iQ; @@ -125,19 +126,19 @@ namespace Grid { 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) + - timesI(toComplex(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0))); + + emiu * (toComplex(16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) + + timesI(toComplex(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0))); - r11 = (toComplex(2.0*unitReal) + timesI(toComplex(4.0*u)))* e2iu - + emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) + - timesI(toComplex(2.0*u*cosw + 6.0*u*xi0))); + r11 = (toComplex(2.0*LatticeUnitReal) + timesI(toComplex(4.0*u)))* e2iu + + emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) + + timesI(toComplex(2.0*u*cosw + 6.0*u*xi0))); - r21 = timesI(toComplex(2.0*unitReal)) * e2iu - + emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0))); + r21 = 2.0*timesI(e2iu) + + emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0))); r02 = -2.0 * e2iu + emiu * (toComplex(-8.0*u2*xi0) + - timesI(toComplex(2.0*u*(cosw + xi0 + 3.0*u2*xi1)))); + timesI(toComplex(2.0*u*(cosw + xi0 + 3.0*u2*xi1)))); r12 = emiu * (toComplex(2.0*u*xi0) + timesI(toComplex(-cosw - xi0 + 3.0*u2*xi1))); @@ -172,34 +173,34 @@ namespace Grid { GaugeLinkField USQ = USigmap * timesMinusI(iQ); GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 + - f1 * USigmap + f2 * QUS + f2 * USQ; + f1 * USigmap + f2 * QUS + f2 * USQ; iLambda = Ta(iGamma); - - } + +} //==================================================================== - public: +public: GaugeField* ThinLinks; /*!< @brief Pointer to the thin links configuration */ /*! @brief Standard constructor */ - SmearedConfiguration(GridCartesian * UGrid, - unsigned int Nsmear, - Smear_Stout& Stout): + SmearedConfiguration(GridCartesian * UGrid, + unsigned int Nsmear, + Smear_Stout& Stout): smearingLevels(Nsmear), - StoutSmearing(Stout), - ThinLinks(NULL){ - for (unsigned int i=0; i< smearingLevels; ++i) - SmearedSet.push_back(*(new GaugeField(UGrid))); + StoutSmearing(Stout), + ThinLinks(NULL){ + for (unsigned int i=0; i< smearingLevels; ++i) + SmearedSet.push_back(*(new GaugeField(UGrid))); } /*! For just thin links */ - SmearedConfiguration(): + SmearedConfiguration(): smearingLevels(0), - StoutSmearing(), - SmearedSet(), - ThinLinks(NULL){} + StoutSmearing(), + SmearedSet(), + ThinLinks(NULL){} // attach the smeared routines to the thin links U and fill the smeared set @@ -207,59 +208,59 @@ namespace Grid { //==================================================================== void smeared_force(GaugeField& SigmaTilde) const{ - if (smearingLevels > 0){ + if (smearingLevels > 0){ GaugeField force = SigmaTilde;//actually = U*SigmaTilde GaugeLinkField tmp_mu(SigmaTilde._grid); - for (int mu = 0; mu < Nd; mu++){ + 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); + 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,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); + tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); + pokeLorentz(SigmaTilde, tmp_mu, mu); } }// if smearingLevels = 0 do nothing - } +} //==================================================================== - - GaugeField& get_SmearedU(){ + +GaugeField& get_SmearedU(){ return SmearedSet[smearingLevels-1]; - } - - GaugeField& get_U(bool smeared=false) { +} + +GaugeField& get_U(bool smeared=false) { // get the config, thin links by default if (smeared){ - 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; - } + 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{ - RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); - std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; - return *ThinLinks;} - } - - }; - - - } - + 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 5b5e4ba3..ff7a1cb8 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -5,64 +5,64 @@ #ifndef STOUT_SMEAR_ #define STOUT_SMEAR_ -namespace Grid { - namespace QCD { + namespace Grid { + namespace QCD { /*! @brief Stout smearing of link variable. */ template - class Smear_Stout: public Smear { - private: - const std::vector d_rho; - const Smear < Gimpl > * SmearBase; - + class Smear_Stout: public Smear { + private: + const std::vector d_rho; + const Smear < Gimpl > * SmearBase; + + + + public: + INHERIT_GIMPL_TYPES(Gimpl) + + Smear_Stout(Smear < Gimpl >* base):SmearBase(base){ + static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); + } - - public: - INHERIT_GIMPL_TYPES(Gimpl) - - Smear_Stout(Smear < Gimpl >* base):SmearBase(base){ - static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); - } - /*! Default constructor */ - Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){ - static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); - } - - ~Smear_Stout(){} - - void smear(GaugeField& u_smr,const GaugeField& U) const{ - GaugeField C(U._grid); - GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid); - - std::cout<< GridLogDebug << "Stout smearing started\n"; - + Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){ + static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); + } + + ~Smear_Stout(){} + + void smear(GaugeField& u_smr,const GaugeField& U) const{ + GaugeField C(U._grid); + GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid); + + std::cout<< GridLogDebug << "Stout smearing started\n"; + //Smear the configurations - SmearBase->smear(C, U); - for (int mu = 0; musmear(C, U); + for (int mu = 0; muderivative(SigmaTerm, iLambda, Gauge); - } - - - void BaseSmear(GaugeField& C, - const GaugeField& U) const{ +} + + +void BaseSmear(GaugeField& C, + const GaugeField& U) const{ SmearBase->smear(C, U); - } - - void exponentiate_iQ(GaugeLinkField& e_iQ, - const GaugeLinkField& iQ) const{ +} + +void exponentiate_iQ(GaugeLinkField& e_iQ, + const GaugeLinkField& iQ) const{ // Put this outside // only valid for SU(3) matrices @@ -87,11 +87,11 @@ namespace Grid { set_uw(u, w, iQ2, iQ3); set_fj(f0, f1, f2, u, w); e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2; - }; +}; - void set_uw(LatticeReal& u, LatticeReal& w, - GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{ +void set_uw(LatticeReal& u, LatticeReal& w, + GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{ Real one_over_three = 1.0/3.0; Real one_over_two = 1.0/2.0; @@ -99,18 +99,20 @@ namespace Grid { LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid); // sign in c0 from the conventions on the Ta - c0 = - toReal(imag(trace(iQ3))) * one_over_three; +// c0 = - toReal(imag(trace(iQ3))) * one_over_three; + c0 = - toReal(real(timesMinusI(trace(iQ3)))) * one_over_three; //temporary, FIX the bug in imag c1 = - toReal(real(trace(iQ2))) * one_over_two; tmp = c1 * one_over_three; - c0max = 2.0 * pow(tmp, 1.5); + c0max = 2.0 * pow(tmp, 1.5); theta = acos(c0/c0max); u = sqrt(tmp) * cos( theta * one_over_three); w = sqrt(c1) * sin( theta * one_over_three); - } - void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2, - const LatticeReal& u, const LatticeReal& w) const{ +} + +void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2, + const LatticeReal& u, const LatticeReal& w) const{ GridBase *grid = u._grid; LatticeReal xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid); @@ -124,41 +126,42 @@ namespace Grid { cosw = cos(w); ixi0 = timesI(toComplex(xi0)); - 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))); h0 = e2iu * toComplex(u2 - w2) + emiu *( toComplex(8.0*u2*cosw) + - toComplex(2.0*u*(3.0*u2 + w2))*ixi0); - + toComplex(2.0*u*(3.0*u2 + w2))*ixi0); + h1 = toComplex(2.0*u) * e2iu - emiu*( toComplex(2.0*u*cosw) - - toComplex(3.0*u2-w2)*ixi0); + toComplex(3.0*u2-w2)*ixi0); h2 = e2iu - emiu * (toComplex(cosw) + toComplex(3.0*u)*ixi0); - tmp = 9.0*u2 - w2; fden = toComplex(pow(tmp, -1.0)); f0 = h0 * fden; f1 = h1 * fden; f2 = h2 * fden; - } + + +} - LatticeReal func_xi0(const LatticeReal& w) const{ +LatticeReal func_xi0(const LatticeReal& w) const{ // Define a function to do the check //if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: "<< w <<"\n"; return sin(w)/w; - } +} - LatticeReal func_xi1(const LatticeReal& w) const{ +LatticeReal func_xi1(const LatticeReal& w) const{ // Define a function to do the check //if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n"; return cos(w)/(w*w) - sin(w)/(w*w*w); - } - - }; +} - } +}; + +} } #endif