portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 05:54:32 +00:00
Code Releases Activity

Debugged set_fj,

Browse Source 
to be fixed: BUG in imag()
This commit is contained in:
Guido Cossu
2016-07-01 16:06:20 +01:00
parent 565e9329ba
commit 092fa0d8da
5 changed files with 233 additions and 228 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
lib/Log.h
View File

@@ -63,7 +63,6 @@ namespace Grid {
colour["WHITE"] ="\033[37m"; colour["WHITE"] ="\033[37m";
colour["NORMAL"] ="\033[0;39m"; colour["NORMAL"] ="\033[0;39m";
} else { } else {
std::cout << "Switching off colours\n";
colour["BLACK"] =""; colour["BLACK"] ="";
colour["RED"] =""; colour["RED"] ="";
colour["GREEN"] =""; colour["GREEN"] ="";

6
lib/qcd/hmc/integrators/Integrator.h
View File

@@ -117,10 +117,10 @@ namespace Grid{
GaugeField force(U._grid); GaugeField force(U._grid);
GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); 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 as[level].actions.at(a)->deriv(Us,force); // deriv should not include Ta
std::cout<<GridLogIntegrator<< "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl; std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl;
if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
force = Ta(force); force = Ta(force);
std::cout<<GridLogIntegrator<< "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl; std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl;
Mom = Mom - force*ep; Mom = Mom - force*ep;
} }
} }
@@ -130,7 +130,7 @@ namespace Grid{
t_U+=ep; t_U+=ep;
int fl = levels-1; int fl = levels-1;
std::cout<<GridLogIntegrator<<" "<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl; std::cout<< GridLogIntegrator <<" "<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
} }
void update_U(GaugeField &Mom, GaugeField&U, double ep){ void update_U(GaugeField &Mom, GaugeField&U, double ep){

116
lib/qcd/smearing/APEsmearing.h
View File

@@ -5,69 +5,71 @@
#ifndef APE_SMEAR_ #ifndef APE_SMEAR_
#define APE_SMEAR_ #define APE_SMEAR_
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
/*! @brief APE type smearing of link variables. */ /*! @brief APE type smearing of link variables. */
template <class Gimpl> template <class Gimpl>
class Smear_APE: public Smear<Gimpl>{ class Smear_APE: public Smear<Gimpl>{
private: private:
const std::vector<double> rho;/*!< Array of weights */ const std::vector<double> rho;/*!< Array of weights */
//This member must be private - we do not want to control from outside //This member must be private - we do not want to control from outside
std::vector<double> set_rho(const double common_rho)const { std::vector<double> set_rho(const double common_rho)const {
std::vector<double> res; std::vector<double> res;
for(int mn=0; mn<Nd*Nd; ++mn) res.push_back(common_rho);
for(int mu=0; mu<Nd; ++mu) res[mu + mu*Nd] = 0.0;
return res;
}
public: for(int mn=0; mn<Nd*Nd; ++mn) res.push_back(common_rho);
for(int mu=0; mu<Nd; ++mu) res[mu + mu*Nd] = 0.0;
return res;
}
public:
// Defines the gauge field types // Defines the gauge field types
INHERIT_GIMPL_TYPES(Gimpl) INHERIT_GIMPL_TYPES(Gimpl)
// Constructors and destructors // Constructors and destructors
Smear_APE(const std::vector<double>& rho_):rho(rho_){} Smear_APE(const std::vector<double>& rho_):rho(rho_){}
Smear_APE(double rho_val):rho(set_rho(rho_val)){} Smear_APE(double rho_val):rho(set_rho(rho_val)){}
Smear_APE():rho(set_rho(1.0)){} Smear_APE():rho(set_rho(1.0)){}
~Smear_APE(){} ~Smear_APE(){}
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////
void smear(GaugeField& u_smr, const GaugeField& U)const{ void smear(GaugeField& u_smr, const GaugeField& U)const{
GridBase *grid = U._grid; GridBase *grid = U._grid;
double d_rho; double d_rho;
GaugeLinkField Cup(grid), tmp_stpl(grid); GaugeLinkField Cup(grid), tmp_stpl(grid);
WilsonLoops<Gimpl> WL; WilsonLoops<Gimpl> WL;
u_smr = zero; u_smr = zero;
for(int mu=0; mu<Nd; ++mu){ for(int mu=0; mu<Nd; ++mu){
Cup = zero; Cup = zero;
for(int nu=0; nu<Nd; ++nu){ for(int nu=0; nu<Nd; ++nu){
if (nu != mu) { if (nu != mu) {
d_rho = rho[mu + Nd * nu]; d_rho = rho[mu + Nd * nu];
// get the staple in direction mu, nu // get the staple in direction mu, nu
WL.Staple(tmp_stpl, U, mu, nu); //nb staple conventions of IroIro and Grid differ by a dagger WL.Staple(tmp_stpl, U, mu, nu); //nb staple conventions of IroIro and Grid differ by a dagger
Cup += tmp_stpl*d_rho; Cup += tmp_stpl*d_rho;
}
} }
}
// save the Cup link-field on the u_smr gauge-field // save the Cup link-field on the u_smr gauge-field
pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag
} }
} }
////////////////////////////////////////////////////////////////////////////////
void derivative(GaugeField& SigmaTerm,
const GaugeField& iLambda,
const GaugeField& U)const{
////////////////////////////////////////////////////////////////////////////////
void derivative(GaugeField& SigmaTerm,
const GaugeField& iLambda,
const GaugeField& U)const{
// Reference // Reference
// Morningstar, Peardon, Phys.Rev.D69,054501(2004) // Morningstar, Peardon, Phys.Rev.D69,054501(2004)
// Equation 75 // Equation 75
// Computing Sigma_mu, derivative of S[fat links] with respect to the thin links
// Output SigmaTerm
GridBase *grid = U._grid; GridBase *grid = U._grid;
int vol = U._grid->gSites(); int vol = U._grid->gSites();
WilsonLoops<Gimpl> WL; WilsonLoops<Gimpl> WL;
@@ -76,27 +78,27 @@ namespace Grid {
GaugeLinkField U_mu(grid), U_nu(grid); GaugeLinkField U_mu(grid), U_nu(grid);
GaugeLinkField sh_field(grid), temp_Sigma(grid); GaugeLinkField sh_field(grid), temp_Sigma(grid);
Real rho_munu, rho_numu; Real rho_munu, rho_numu;
for(int mu = 0; mu < Nd; ++mu){
U_mu = PeekIndex<LorentzIndex>( U, mu);
iLambda_mu = PeekIndex<LorentzIndex>(iLambda, mu);
for(int nu = 0; nu < Nd; ++nu){
if(nu==mu) continue;
U_nu = PeekIndex<LorentzIndex>( U, nu);
iLambda_nu = PeekIndex<LorentzIndex>(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<LorentzIndex>( U, mu);
iLambda_mu = PeekIndex<LorentzIndex>(iLambda, mu);
for(int nu = 0; nu < Nd; ++nu){
if(nu==mu) continue;
U_nu = PeekIndex<LorentzIndex>( U, nu);
iLambda_nu = PeekIndex<LorentzIndex>(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) //-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? sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
temp_Sigma = rho_numu*sh_field*staple; temp_Sigma = rho_numu*sh_field*staple;
//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x) //r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu); Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
@@ -119,10 +121,10 @@ namespace Grid {
sh_field = Cshift(temp_Sigma, nu, -1); sh_field = Cshift(temp_Sigma, nu, -1);
Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu); Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
}
} }
} }
}; }
};

203
lib/qcd/smearing/GaugeConfiguration.h
View File

@@ -6,10 +6,10 @@
#ifndef GAUGE_CONFIG_ #ifndef GAUGE_CONFIG_
#define GAUGE_CONFIG_ #define GAUGE_CONFIG_
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
/*! /*!
@brief Smeared configuration container @brief Smeared configuration container
@@ -23,42 +23,43 @@ namespace Grid {
*/ */
template <class Gimpl> template <class Gimpl>
class SmearedConfiguration { class SmearedConfiguration {
public: public:
INHERIT_GIMPL_TYPES(Gimpl) ; INHERIT_GIMPL_TYPES(Gimpl) ;
private: private:
const unsigned int smearingLevels; const unsigned int smearingLevels;
Smear_Stout<Gimpl> StoutSmearing; Smear_Stout<Gimpl> StoutSmearing;
std::vector<GaugeField> SmearedSet; std::vector<GaugeField> SmearedSet;
// Member functions // Member functions
//==================================================================== //====================================================================
void fill_smearedSet(GaugeField& U){ void fill_smearedSet(GaugeField& U){
ThinLinks = &U; //attach the smearing routine to the field U ThinLinks = &U; //attach the smearing routine to the field U
//check the pointer is not null //check the pointer is not null
if (ThinLinks==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){ if (smearingLevels > 0){
std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
GaugeField previous_u(ThinLinks->_grid); 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];
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(previous_u); previous_u = *ThinLinks;
std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl; for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
previous_u = SmearedSet[smearLvl];
} // For debug purposes
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(previous_u);
std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl;
}
} }
} }
//==================================================================== //====================================================================
GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
const GaugeField& GaugeK) const{ const GaugeField& GaugeK) const{
GridBase *grid = GaugeK._grid; GridBase *grid = GaugeK._grid;
GaugeField C(grid), SigmaK(grid), iLambda(grid); GaugeField C(grid), SigmaK(grid), iLambda(grid);
GaugeLinkField iLambda_mu(grid); GaugeLinkField iLambda_mu(grid);
@@ -69,27 +70,27 @@ namespace Grid {
StoutSmearing.BaseSmear(C, GaugeK); StoutSmearing.BaseSmear(C, GaugeK);
for (int mu = 0; mu < Nd; mu++){ for (int mu = 0; mu < Nd; mu++){
Cmu = peekLorentz( C,mu); Cmu = peekLorentz( C,mu);
GaugeKmu = peekLorentz(GaugeK,mu); GaugeKmu = peekLorentz(GaugeK,mu);
SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu); SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu);
iQ = Ta(Cmu*adj(GaugeKmu)); iQ = Ta(Cmu*adj(GaugeKmu));
set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu); set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu);
pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu); pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu);
pokeLorentz(iLambda, iLambda_mu, mu); pokeLorentz(iLambda, iLambda_mu, mu);
} }
StoutSmearing.derivative(SigmaK, iLambda, GaugeK);// derivative of SmearBase StoutSmearing.derivative(SigmaK, iLambda, GaugeK);// derivative of SmearBase
return SigmaK; return SigmaK;
} }
/*! @brief Returns smeared configuration at level 'Level' */ /*! @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]; return SmearedSet[Level];
} }
//==================================================================== //====================================================================
void set_iLambda(GaugeLinkField& iLambda, void set_iLambda(GaugeLinkField& iLambda,
GaugeLinkField& e_iQ, GaugeLinkField& e_iQ,
const GaugeLinkField& iQ, const GaugeLinkField& iQ,
const GaugeLinkField& Sigmap, const GaugeLinkField& Sigmap,
const GaugeLinkField& GaugeK)const{ const GaugeLinkField& GaugeK)const{
GridBase *grid = iQ._grid; GridBase *grid = iQ._grid;
GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid); GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid);
GaugeLinkField unity(grid); GaugeLinkField unity(grid);
@@ -103,9 +104,9 @@ namespace Grid {
LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid); LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid);
LatticeComplex r22(grid), tr1(grid), tr2(grid); LatticeComplex r22(grid), tr1(grid), tr2(grid);
LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(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 // Exponential
iQ2 = iQ * iQ; iQ2 = iQ * iQ;
@@ -125,19 +126,19 @@ namespace Grid {
e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(sin(2.0*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 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) + + 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))); timesI(toComplex(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0)));
r11 = (toComplex(2.0*unitReal) + timesI(toComplex(4.0*u)))* e2iu r11 = (toComplex(2.0*LatticeUnitReal) + timesI(toComplex(4.0*u)))* e2iu
+ emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) + + emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) +
timesI(toComplex(2.0*u*cosw + 6.0*u*xi0))); timesI(toComplex(2.0*u*cosw + 6.0*u*xi0)));
r21 = timesI(toComplex(2.0*unitReal)) * e2iu r21 = 2.0*timesI(e2iu)
+ emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0))); + emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0)));
r02 = -2.0 * e2iu + emiu * (toComplex(-8.0*u2*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))); 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 USQ = USigmap * timesMinusI(iQ);
GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 + GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 +
f1 * USigmap + f2 * QUS + f2 * USQ; f1 * USigmap + f2 * QUS + f2 * USQ;
iLambda = Ta(iGamma); iLambda = Ta(iGamma);
} }
//==================================================================== //====================================================================
public: public:
GaugeField* ThinLinks; /*!< @brief Pointer to the thin GaugeField* ThinLinks; /*!< @brief Pointer to the thin
links configuration */ links configuration */
/*! @brief Standard constructor */ /*! @brief Standard constructor */
SmearedConfiguration(GridCartesian * UGrid, SmearedConfiguration(GridCartesian * UGrid,
unsigned int Nsmear, unsigned int Nsmear,
Smear_Stout<Gimpl>& Stout): Smear_Stout<Gimpl>& Stout):
smearingLevels(Nsmear), smearingLevels(Nsmear),
StoutSmearing(Stout), StoutSmearing(Stout),
ThinLinks(NULL){ ThinLinks(NULL){
for (unsigned int i=0; i< smearingLevels; ++i) for (unsigned int i=0; i< smearingLevels; ++i)
SmearedSet.push_back(*(new GaugeField(UGrid))); SmearedSet.push_back(*(new GaugeField(UGrid)));
} }
/*! For just thin links */ /*! For just thin links */
SmearedConfiguration(): SmearedConfiguration():
smearingLevels(0), smearingLevels(0),
StoutSmearing(), StoutSmearing(),
SmearedSet(), SmearedSet(),
ThinLinks(NULL){} ThinLinks(NULL){}
// attach the smeared routines to the thin links U and fill the smeared set // 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{ void smeared_force(GaugeField& SigmaTilde) const{
if (smearingLevels > 0){ if (smearingLevels > 0){
GaugeField force = SigmaTilde;//actually = U*SigmaTilde GaugeField force = SigmaTilde;//actually = U*SigmaTilde
GaugeLinkField tmp_mu(SigmaTilde._grid); GaugeLinkField tmp_mu(SigmaTilde._grid);
for (int mu = 0; mu < Nd; mu++){ for (int mu = 0; mu < Nd; mu++){
// to get SigmaTilde // to get SigmaTilde
tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu); tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
pokeLorentz(force, tmp_mu, mu); pokeLorentz(force, tmp_mu, mu);
} }
for(int ismr = smearingLevels - 1; ismr > 0; --ismr) 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); force = AnalyticSmearedForce(force,*ThinLinks);
for (int mu = 0; mu < Nd; mu++){ for (int mu = 0; mu < Nd; mu++){
tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
pokeLorentz(SigmaTilde, tmp_mu, mu); pokeLorentz(SigmaTilde, tmp_mu, mu);
} }
}// if smearingLevels = 0 do nothing }// if smearingLevels = 0 do nothing
} }
//==================================================================== //====================================================================
GaugeField& get_SmearedU(){ GaugeField& get_SmearedU(){
return SmearedSet[smearingLevels-1]; return SmearedSet[smearingLevels-1];
} }
GaugeField& get_U(bool smeared=false) { GaugeField& get_U(bool smeared=false) {
// get the config, thin links by default // get the config, thin links by default
if (smeared){ if (smeared){
if (smearingLevels){ if (smearingLevels){
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(SmearedSet[smearingLevels-1]); RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(SmearedSet[smearingLevels-1]);
std::cout<< GridLogDebug << "getting U Plaq: " << impl_plaq<< std::endl; std::cout<< GridLogDebug << "getting U Plaq: " << impl_plaq<< std::endl;
return get_SmearedU(); return get_SmearedU();
} }
else { else {
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks); RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks);
std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl;
return *ThinLinks; return *ThinLinks;
} }
} }
else{ else{
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks); RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks);
std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl;
return *ThinLinks;} return *ThinLinks;}
} }
}; };
} }
} }

135
lib/qcd/smearing/StoutSmearing.h
View File

@@ -5,64 +5,64 @@
#ifndef STOUT_SMEAR_ #ifndef STOUT_SMEAR_
#define STOUT_SMEAR_ #define STOUT_SMEAR_
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
/*! @brief Stout smearing of link variable. */ /*! @brief Stout smearing of link variable. */
template <class Gimpl> template <class Gimpl>
class Smear_Stout: public Smear<Gimpl> { class Smear_Stout: public Smear<Gimpl> {
private: private:
const std::vector<double> d_rho; const std::vector<double> d_rho;
const Smear < Gimpl > * SmearBase; 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 */ /*! Default constructor */
Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){ Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){
static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
} }
~Smear_Stout(){} ~Smear_Stout(){}
void smear(GaugeField& u_smr,const GaugeField& U) const{ void smear(GaugeField& u_smr,const GaugeField& U) const{
GaugeField C(U._grid); GaugeField C(U._grid);
GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid); GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid);
std::cout<< GridLogDebug << "Stout smearing started\n"; std::cout<< GridLogDebug << "Stout smearing started\n";
//Smear the configurations //Smear the configurations
SmearBase->smear(C, U); SmearBase->smear(C, U);
for (int mu = 0; mu<Nd; mu++){ for (int mu = 0; mu<Nd; mu++){
tmp = peekLorentz(C,mu); tmp = peekLorentz(C,mu);
Umu = peekLorentz(U,mu); Umu = peekLorentz(U,mu);
iq_mu = Ta(tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper iq_mu = Ta(tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
exponentiate_iQ(tmp, iq_mu); exponentiate_iQ(tmp, iq_mu);
pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu
} }
std::cout<< GridLogDebug << "Stout smearing completed\n"; std::cout<< GridLogDebug << "Stout smearing completed\n";
} }
void derivative(GaugeField& SigmaTerm, void derivative(GaugeField& SigmaTerm,
const GaugeField& iLambda, const GaugeField& iLambda,
const GaugeField& Gauge) const{ const GaugeField& Gauge) const{
SmearBase->derivative(SigmaTerm, iLambda, Gauge); SmearBase->derivative(SigmaTerm, iLambda, Gauge);
} }
void BaseSmear(GaugeField& C, void BaseSmear(GaugeField& C,
const GaugeField& U) const{ const GaugeField& U) const{
SmearBase->smear(C, U); SmearBase->smear(C, U);
} }
void exponentiate_iQ(GaugeLinkField& e_iQ, void exponentiate_iQ(GaugeLinkField& e_iQ,
const GaugeLinkField& iQ) const{ const GaugeLinkField& iQ) const{
// Put this outside // Put this outside
// only valid for SU(3) matrices // only valid for SU(3) matrices
@@ -87,11 +87,11 @@ namespace Grid {
set_uw(u, w, iQ2, iQ3); set_uw(u, w, iQ2, iQ3);
set_fj(f0, f1, f2, u, w); set_fj(f0, f1, f2, u, w);
e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2; e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
}; };
void set_uw(LatticeReal& u, LatticeReal& w, void set_uw(LatticeReal& u, LatticeReal& w,
GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{ GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
Real one_over_three = 1.0/3.0; Real one_over_three = 1.0/3.0;
Real one_over_two = 1.0/2.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); LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
// sign in c0 from the conventions on the Ta // 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; c1 = - toReal(real(trace(iQ2))) * one_over_two;
tmp = c1 * one_over_three; tmp = c1 * one_over_three;
c0max = 2.0 * pow(tmp, 1.5); c0max = 2.0 * pow(tmp, 1.5);
theta = acos(c0/c0max); theta = acos(c0/c0max);
u = sqrt(tmp) * cos( theta * one_over_three); u = sqrt(tmp) * cos( theta * one_over_three);
w = sqrt(c1) * sin( 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; GridBase *grid = u._grid;
LatticeReal xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid); LatticeReal xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid);
@@ -124,41 +126,42 @@ namespace Grid {
cosw = cos(w); cosw = cos(w);
ixi0 = timesI(toComplex(xi0)); ixi0 = timesI(toComplex(xi0));
emiu = toComplex(cos(u)) - timesI(toComplex(u)); emiu = toComplex(cos(u)) - timesI(toComplex(sin(u)));
e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u)); e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(sin(2.0*u)));
h0 = e2iu * toComplex(u2 - w2) + emiu *( toComplex(8.0*u2*cosw) + 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) - 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); h2 = e2iu - emiu * (toComplex(cosw) + toComplex(3.0*u)*ixi0);
tmp = 9.0*u2 - w2; tmp = 9.0*u2 - w2;
fden = toComplex(pow(tmp, -1.0)); fden = toComplex(pow(tmp, -1.0));
f0 = h0 * fden; f0 = h0 * fden;
f1 = h1 * fden; f1 = h1 * fden;
f2 = h2 * 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 // Define a function to do the check
//if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: "<< w <<"\n"; //if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: "<< w <<"\n";
return sin(w)/w; 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 // Define a function to do the check
//if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n"; //if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n";
return cos(w)/(w*w) - sin(w)/(w*w*w); return cos(w)/(w*w) - sin(w)/(w*w*w);
} }
};
} };
}
} }
#endif #endif