From c1b1b89d17e8fb06bb81e61691356071bf2afe77 Mon Sep 17 00:00:00 2001 From: neo Date: Fri, 19 Feb 2016 17:15:27 +0900 Subject: [PATCH] More on smearing routines, writing APEsmear (dev) --- lib/qcd/smearing/APEsmearing.h | 150 ++++++++++++++++++++++++++ lib/qcd/smearing/BaseSmearing.h | 17 +++ lib/qcd/smearing/GaugeConfiguration.h | 16 ++- lib/qcd/smearing/StoutSmearing.h | 32 ++++++ lib/qcd/utils/WilsonLoops.h | 83 ++++++++++++++ 5 files changed, 294 insertions(+), 4 deletions(-) create mode 100644 lib/qcd/smearing/APEsmearing.h create mode 100644 lib/qcd/smearing/BaseSmearing.h create mode 100644 lib/qcd/smearing/StoutSmearing.h diff --git a/lib/qcd/smearing/APEsmearing.h b/lib/qcd/smearing/APEsmearing.h new file mode 100644 index 00000000..e2a52b56 --- /dev/null +++ b/lib/qcd/smearing/APEsmearing.h @@ -0,0 +1,150 @@ +/*! + @brief Declaration of Smear_APE class for APE smearing +*/ + +#ifndef APE_SMEAR_ +#define APE_SMEAR_ + +/*! @brief APE type smearing of link variables. */ + +template +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)const { + std::vector res; + + 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(){} + + void smear(GaugeField& u_smr, const GaugeField& U)const{ + double d_rho; + GaugeLinkField Cup, tmp_stpl; + WilsonLoops WL; + u_smr = zero; + + for(int mu=0; mu WL; + GaugeLinkField staple, u_tmp, iLambda_mu, iLambda_nu; + GaugeLinkField U_mu, U_nu; + GaugeLinkField sh_field ; + GaugeLinkField temp_Sigma; + + SU::Matrix temp_mat, temp_mat2; + Real rho_munu, rho_numu; + + // to be completed + int Nvol = CommonPrms::instance()->Nvol(); + + 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 = adj(staple)*iLambda_nu; + temp_Sigma *= - rho_numu; + //-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x) + SigmaTerm ................. + + for (int site = 0; site < Nvol; ++site){ + temp_mat = mat_dag(staple,site) * mat(iLambda_nu,site); + temp_mat *= - rho_numu; + AddMat(SigmaTerm, temp_mat, site, mu); + } + sh_field = shiftField(iLambda_nu, mu, Forward()); + + for (int site = 0; site < Nvol; ++site){ + temp_mat = mat(sh_field,site) * mat_dag(staple,site); + temp_mat *= rho_numu; + AddMat(SigmaTerm, temp_mat, site, mu); + }//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x) + + sh_field = shiftField(iLambda_mu, nu, Forward()); + + for (int site = 0; site < Nvol; ++site){ + temp_mat = mat(U_nu,site) * mat(sh_field,site) * mat_dag(U_nu,site); + temp_mat = mat_dag(staple,site) * temp_mat; + temp_mat *= - rho_munu; + AddMat(SigmaTerm, temp_mat, site, mu); + }//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x) + + staple = 0.0; + sh_field = shiftField(U_nu, mu, Forward()); + + for (int site = 0; site < Nvol; ++site){ + temp_mat2 = mat_dag(sh_field,site) * mat_dag(U_mu,site); + temp_mat = temp_mat2 * mat(iLambda_mu,site) * mat(U_nu,site); + temp_mat *= - rho_munu; + AddMat(staple, temp_mat, site); + temp_mat = temp_mat2 * mat(iLambda_nu,site) * mat(U_nu,site); + temp_mat *= rho_numu; + AddMat(staple, temp_mat, site); + } + + for (int site = 0; site < Nvol; ++site){ + temp_mat = mat_dag(U_nu,site) * mat(iLambda_nu,site); + SetMat(u_tmp, temp_mat, site); + } + + sh_field = shiftField(u_tmp, mu, Forward()); + + for (int site = 0; site < Nvol; ++site){ + temp_mat = mat(sh_field,site) * mat_dag(U_mu,site) * mat(U_nu,site); + temp_mat *= - rho_numu; + AddMat(staple, temp_mat, site); + } + + sh_field = shiftField(staple, nu, Backward()); + + AddSlice(SigmaTerm, sh_field, mu); + } + } + + */ + } + + + + +}; + +#endif diff --git a/lib/qcd/smearing/BaseSmearing.h b/lib/qcd/smearing/BaseSmearing.h new file mode 100644 index 00000000..377bcbc9 --- /dev/null +++ b/lib/qcd/smearing/BaseSmearing.h @@ -0,0 +1,17 @@ +/* + @brief Declares base smearing class Smear + */ +#ifndef BASE_SMEAR_ +#define BASE_SMEAR_ + +template +class Smear{ +public: + INHERIT_GIMPL_TYPES(Gimpl) + + virtual ~Smear(){} + virtual void smear (GaugeField&,const GaugeField&)const = 0; + virtual void derivative(GaugeField&, + const GaugeField&,const GaugeField&) const = 0; +}; +#endif diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h index d040b739..5da5d8fa 100644 --- a/lib/qcd/smearing/GaugeConfiguration.h +++ b/lib/qcd/smearing/GaugeConfiguration.h @@ -77,11 +77,15 @@ namespace Grid { void set_GaugeField(){ fill_smearedSet(); } void smeared_force(GaugeField&) const; - GaugeField& get_current_conf() const; - GaugeField& select_conf(bool smeared) const { + GaugeField* get_SmearedU() const{ + return const_cast(&(SmearedSet[smearingLevels-1])); + } + + GaugeField* get_U(bool smeared=false) const { + // get the config, thin links by default if (smeared){ - if (smearingLevels) return get_current_conf(); - else return ThinLinks; + if (smearingLevels) return get_SmearedU(); + else return ThinLinks; } else return ThinLinks; } @@ -93,4 +97,8 @@ namespace Grid { } + + + + #endif diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h new file mode 100644 index 00000000..4b372d64 --- /dev/null +++ b/lib/qcd/smearing/StoutSmearing.h @@ -0,0 +1,32 @@ +/* + @file stoutSmear.hpp + @brief Declares Stout smearing class +*/ +#ifndef STOUT_SMEAR_ +#define STOUT_SMEAR_ + +/*! @brief Stout smearing of link variable. */ +template +class Smear_Stout: public Smear { +private: + const std::valarray d_rho; + const Smear* SmearBase; + + double func_xi0(double w) const; +public: + INHERIT_GIMPL_TYPES(Gimpl) + Smear_Stout(Smear* base):SmearBase(base){} + + /*! Default constructor */ + Smear_Stout():SmearBase(new Smear_APE()){} + + ~Smear_Stout(){} + + void smear(GaugeField&,const GaugeField&) const; + void BaseSmear(GaugeField&, const GaugeField&) const; + void derivative(GaugeField&, const GaugeField&, const GaugeField&) const; + void exponentiate_iQ(GaugeField&, const GaugeField&) const; + +}; + +#endif diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index bab6dec9..8272f1dc 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -171,6 +171,89 @@ public: } } + ////////////////////////////////////////////////// + // the sum over all staples on each site in direction mu,nu + ////////////////////////////////////////////////// + static void Staple(GaugeMat &staple,const GaugeLorentz &Umu,int mu, in nu){ + + GridBase *grid = Umu._grid; + + std::vector U(4,grid); + for(int d=0;d(Umu,d); + } + staple = zero; + GaugeMat tmp(grid); + + + if(nu != mu) { + + // mu + // ^ + // |__> nu + + // __ + // | + // __| + // + + staple+=Gimpl::ShiftStaple( + Gimpl::CovShiftForward (U[nu],nu, + Gimpl::CovShiftBackward(U[mu],mu, + Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu); + + // __ + // | + // |__ + // + // + staple+=Gimpl::ShiftStaple( + Gimpl::CovShiftBackward(U[nu],nu, + Gimpl::CovShiftBackward(U[mu],mu,U[nu])),mu); + + } + } + + + ////////////////////////////////////////////////// + // the sum over all staples on each site in direction mu,nu, upper part + ////////////////////////////////////////////////// + static void StapleUpper(GaugeMat &staple,const GaugeLorentz &Umu,int mu, in nu){ + + GridBase *grid = Umu._grid; + + std::vector U(4,grid); + for(int d=0;d(Umu,d); + } + staple = zero; + GaugeMat tmp(grid); + + + if(nu != mu) { + + // mu + // ^ + // |__> nu + + // __ + // | + // __| + // + + staple+=Gimpl::ShiftStaple( + Gimpl::CovShiftForward (U[nu],nu, + Gimpl::CovShiftBackward(U[mu],mu, + Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu); + + } + } + + + + + + ////////////////////////////////////////////////////// // Similar to above for rectangle is required //////////////////////////////////////////////////////