portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

Debugged smearing for EOWilson, accepts

Browse Source
This commit is contained in:
Guido Cossu 2016-07-04 15:35:37 +01:00
parent 8dd099267d
commit 0fa66e8f3c
4 changed files with 305 additions and 400 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -44,40 +44,40 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <memory> #include <memory>
namespace Grid{ namespace Grid{
namespace QCD{ namespace QCD{
struct IntegratorParameters{ struct IntegratorParameters{
int Nexp; int Nexp;
int MDsteps; // number of outer steps int MDsteps; // number of outer steps
RealD trajL; // trajectory length RealD trajL; // trajectory length
RealD stepsize; RealD stepsize;
IntegratorParameters(int MDsteps_, IntegratorParameters(int MDsteps_,
RealD trajL_=1.0, RealD trajL_=1.0,
int Nexp_=12): int Nexp_=12):
Nexp(Nexp_), Nexp(Nexp_),
MDsteps(MDsteps_), MDsteps(MDsteps_),
trajL(trajL_), trajL(trajL_),
stepsize(trajL/MDsteps) stepsize(trajL/MDsteps)
{ {
// empty body constructor // empty body constructor
}; };
}; };
/*! @brief Class for Molecular Dynamics management */ /*! @brief Class for Molecular Dynamics management */
template<class GaugeField, class SmearingPolicy> template<class GaugeField, class SmearingPolicy>
class Integrator { class Integrator {
protected: protected:
typedef IntegratorParameters ParameterType; typedef IntegratorParameters ParameterType;
IntegratorParameters Params; IntegratorParameters Params;
const ActionSet<GaugeField> as; const ActionSet<GaugeField> as;
int levels; // int levels; //
double t_U; // Track time passing on each level and for U and for P double t_U; // Track time passing on each level and for U and for P
@ -90,14 +90,14 @@ namespace Grid{
// Should match any legal (SU(n)) gauge field // Should match any legal (SU(n)) gauge field
// Need to use this template to match Ncol to pass to SU<N> class // Need to use this template to match Ncol to pass to SU<N> class
template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){ template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){
typedef Lattice< iScalar< iScalar< iMatrix<vec,Ncol> > > > GaugeLinkField; typedef Lattice< iScalar< iScalar< iMatrix<vec,Ncol> > > > GaugeLinkField;
GaugeLinkField Pmu(P._grid); GaugeLinkField Pmu(P._grid);
Pmu = zero; Pmu = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu); SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
PokeIndex<LorentzIndex>(P, Pmu, mu); PokeIndex<LorentzIndex>(P, Pmu, mu);
}
} }
}
//ObserverList observers; // not yet //ObserverList observers; // not yet
@ -105,126 +105,128 @@ namespace Grid{
// void register_observers(); // void register_observers();
// void notify_observers(); // void notify_observers();
void update_P(GaugeField&U, int level,double ep){ void update_P(GaugeField&U, int level, double ep){
t_P[level]+=ep; t_P[level]+=ep;
update_P(P,U,level,ep); update_P(P,U,level,ep);
std::cout<<GridLogIntegrator<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl; std::cout<<GridLogIntegrator<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
} }
void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){ void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){
for(int a=0; a<as[level].actions.size(); ++a){ // input U actually not used...
GaugeField force(U._grid); for(int a=0; a<as[level].actions.size(); ++a){
GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); GaugeField force(U._grid);
as[level].actions.at(a)->deriv(Us,force); // deriv should not include Ta GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl; as[level].actions.at(a)->deriv(Us,force); // deriv should NOT include Ta
if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
force = Ta(force); std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl;
std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl; if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
Mom = Mom - force*ep; force = Ta(force);
std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl;
Mom -= force*ep;
}
} }
}
void update_U(GaugeField&U, double ep){ void update_U(GaugeField&U, double ep){
update_U(P,U,ep); update_U(P,U,ep);
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){
//rewrite exponential to deal automatically with the lorentz index? //rewrite exponential to deal automatically with the lorentz index?
// GaugeLinkField Umu(U._grid); // GaugeLinkField Umu(U._grid);
// GaugeLinkField Pmu(U._grid); // GaugeLinkField Pmu(U._grid);
for (int mu = 0; mu < Nd; mu++){ for (int mu = 0; mu < Nd; mu++){
auto Umu=PeekIndex<LorentzIndex>(U, mu); auto Umu=PeekIndex<LorentzIndex>(U, mu);
auto Pmu=PeekIndex<LorentzIndex>(Mom, mu); auto Pmu=PeekIndex<LorentzIndex>(Mom, mu);
Umu = expMat(Pmu, ep, Params.Nexp)*Umu; Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
ProjectOnGroup(Umu); ProjectOnGroup(Umu);
PokeIndex<LorentzIndex>(U, Umu, mu); PokeIndex<LorentzIndex>(U, Umu, mu);
} }
// Update the smeared fields, can be implemented as observer // Update the smeared fields, can be implemented as observer
Smearer.set_GaugeField(U); Smearer.set_GaugeField(U);
} }
virtual void step (GaugeField& U,int level, int first,int last)=0;
public: virtual void step (GaugeField& U,int level, int first,int last)=0;
Integrator(GridBase* grid, public:
IntegratorParameters Par,
ActionSet<GaugeField> & Aset, Integrator(GridBase* grid,
SmearingPolicy &Sm): IntegratorParameters Par,
Params(Par), ActionSet<GaugeField> & Aset,
as(Aset), SmearingPolicy &Sm):
P(grid), Params(Par),
levels(Aset.size()), as(Aset),
Smearer(Sm) P(grid),
{ levels(Aset.size()),
t_P.resize(levels,0.0); Smearer(Sm)
t_U=0.0; {
t_P.resize(levels,0.0);
t_U=0.0;
// initialization of smearer delegated outside of Integrator // initialization of smearer delegated outside of Integrator
}; };
virtual ~Integrator(){} virtual ~Integrator(){}
//Initialization of momenta and actions //Initialization of momenta and actions
void refresh(GaugeField& U,GridParallelRNG &pRNG){ void refresh(GaugeField& U,GridParallelRNG &pRNG){
std::cout<<GridLogIntegrator<< "Integrator refresh\n"; std::cout<<GridLogIntegrator<< "Integrator refresh\n";
generate_momenta(P,pRNG); generate_momenta(P,pRNG);
for(int level=0; level< as.size(); ++level){ for(int level=0; level< as.size(); ++level){
for(int actionID=0; actionID<as[level].actions.size(); ++actionID){ for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
// get gauge field from the SmearingPolicy and // get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID // based on the boolean is_smeared in actionID
GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared); GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
as[level].actions.at(actionID)->refresh(Us, pRNG); as[level].actions.at(actionID)->refresh(Us, pRNG);
} }
}
} }
}
// Calculate action // Calculate action
RealD S(GaugeField& U){ RealD S(GaugeField& U){// here also U not used
LatticeComplex Hloc(U._grid); Hloc = zero; LatticeComplex Hloc(U._grid); Hloc = zero;
// Momenta // Momenta
for (int mu=0; mu <Nd; mu++){ for (int mu=0; mu <Nd; mu++){
auto Pmu = PeekIndex<LorentzIndex>(P, mu); auto Pmu = PeekIndex<LorentzIndex>(P, mu);
Hloc -= trace(Pmu*Pmu); Hloc -= trace(Pmu*Pmu);
} }
Complex Hsum = sum(Hloc); Complex Hsum = sum(Hloc);
RealD H = Hsum.real(); RealD H = Hsum.real();
RealD Hterm; RealD Hterm;
std::cout<<GridLogMessage << "Momentum action H_p = "<< H << "\n"; std::cout<<GridLogMessage << "Momentum action H_p = "<< H << "\n";
// Actions // Actions
for(int level=0; level<as.size(); ++level){ for(int level=0; level<as.size(); ++level){
for(int actionID=0; actionID<as[level].actions.size(); ++actionID){ for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
// get gauge field from the SmearingPolicy and // get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID // based on the boolean is_smeared in actionID
GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared); GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
Hterm = as[level].actions.at(actionID)->S(Us); Hterm = as[level].actions.at(actionID)->S(Us);
std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl; std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
H += Hterm; H += Hterm;
} }
} }
return H;
}
void integrate(GaugeField& U){ return H;
}
void integrate(GaugeField& U){
// reset the clocks // reset the clocks
t_U=0; t_U=0;
for(int level=0; level<as.size(); ++level){ for(int level=0; level<as.size(); ++level){
t_P[level]=0; t_P[level]=0;
} }
for(int step=0; step< Params.MDsteps; ++step){ // MD step for(int step=0; step< Params.MDsteps; ++step){ // MD step
int first_step = (step==0); int first_step = (step==0);
int last_step = (step==Params.MDsteps-1); int last_step = (step==Params.MDsteps-1);
this->step(U,0,first_step,last_step); this->step(U,0,first_step,last_step);
} }
// Check the clocks all match on all levels // Check the clocks all match on all levels
@ -237,9 +239,9 @@ namespace Grid{
assert(fabs(t_U-Params.trajL) < 1.0e-6); assert(fabs(t_U-Params.trajL) < 1.0e-6);
} }
}; };
} }
} }
#endif//INTEGRATOR_INCLUDED #endif//INTEGRATOR_INCLUDED

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

@ -13,10 +13,10 @@
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 mn=0; mn<Nd*Nd; ++mn) res.push_back(common_rho);
@ -30,7 +30,7 @@
// Constructors and destructors // Constructors and destructors
Smear_APE(const std::vector<double>& rho_):rho(rho_){} Smear_APE(const std::vector<double>& rho_):rho(rho_){} // check vector size
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(){}
@ -38,7 +38,6 @@
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////
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;
GaugeLinkField Cup(grid), tmp_stpl(grid); GaugeLinkField Cup(grid), tmp_stpl(grid);
WilsonLoops<Gimpl> WL; WilsonLoops<Gimpl> WL;
u_smr = zero; u_smr = zero;
@ -47,21 +46,20 @@
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]; // 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*rho[mu + Nd * nu];
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 see conventions for Staple
pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag }
} }
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void derivative(GaugeField& SigmaTerm, void derivative(GaugeField& SigmaTerm,
const GaugeField& iLambda, const GaugeField& iLambda,
const GaugeField& U)const{ const GaugeField& U)const{
// Reference // Reference
// Morningstar, Peardon, Phys.Rev.D69,054501(2004) // Morningstar, Peardon, Phys.Rev.D69,054501(2004)
@ -69,62 +67,61 @@ void derivative(GaugeField& SigmaTerm,
// Computing Sigma_mu, derivative of S[fat links] with respect to the thin links // Computing Sigma_mu, derivative of S[fat links] with respect to the thin links
// Output SigmaTerm // Output SigmaTerm
GridBase *grid = U._grid; GridBase *grid = U._grid;
int vol = U._grid->gSites();
WilsonLoops<Gimpl> WL;
GaugeLinkField staple(grid), u_tmp(grid);
GaugeLinkField iLambda_mu(grid), iLambda_nu(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){ WilsonLoops<Gimpl> WL;
U_mu = PeekIndex<LorentzIndex>( U, mu); GaugeLinkField staple(grid), u_tmp(grid);
iLambda_mu = PeekIndex<LorentzIndex>(iLambda, mu); GaugeLinkField iLambda_mu(grid), iLambda_nu(grid);
GaugeLinkField U_mu(grid), U_nu(grid);
GaugeLinkField sh_field(grid), temp_Sigma(grid);
Real rho_munu, rho_numu;
for(int nu = 0; nu < Nd; ++nu){ for(int mu = 0; mu < Nd; ++mu){
if(nu==mu) continue; U_mu = peekLorentz( U, mu);
U_nu = PeekIndex<LorentzIndex>( U, nu); iLambda_mu = peekLorentz(iLambda, mu);
iLambda_nu = PeekIndex<LorentzIndex>(iLambda, nu);
rho_munu = rho[mu + Nd * nu]; for(int nu = 0; nu < Nd; ++nu){
rho_numu = rho[nu + Nd * mu]; if(nu==mu) continue;
U_nu = peekLorentz( U, nu);
iLambda_nu = peekLorentz(iLambda, nu);
WL.StapleUpper(staple, U, mu, nu); rho_munu = rho[mu + Nd * nu];
rho_numu = rho[nu + Nd * mu];
temp_Sigma = -rho_numu*staple*iLambda_nu; WL.StapleUpper(staple, U, mu, nu);
//-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity? temp_Sigma = -rho_numu*staple*iLambda_nu; //ok
//-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
temp_Sigma = rho_numu*sh_field*staple; sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
sh_field = Cshift(iLambda_mu, nu, 1); temp_Sigma = rho_numu*sh_field*staple; //ok
//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu); sh_field = Cshift(iLambda_mu, nu, 1);
//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
staple = zero; temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu); //ok
sh_field = Cshift(U_nu, mu, 1); //-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu; staple = zero;
temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu; sh_field = Cshift(U_nu, mu, 1);
u_tmp = adj(U_nu)*iLambda_nu; temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;
sh_field = Cshift(u_tmp, mu, 1); temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu;
temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
sh_field = Cshift(temp_Sigma, nu, -1); u_tmp = adj(U_nu)*iLambda_nu;
Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu); sh_field = Cshift(u_tmp, mu, 1);
temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
} sh_field = Cshift(temp_Sigma, nu, -1);
} Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
}
}; }
}
}
};

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

@ -17,7 +17,7 @@
the HMC update and integrators. the HMC update and integrators.
An "advanced configuration" object that can provide not only the An "advanced configuration" object that can provide not only the
data to store the gauge configuration but also operations to manipulate data to store the gauge configuration but also operations to manipulate
it like smearing. it, like smearing.
It stores a list of smeared configurations. It stores a list of smeared configurations.
*/ */
@ -57,7 +57,7 @@
} }
} }
//==================================================================== //====================================================================
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;
@ -68,12 +68,14 @@ GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
GaugeLinkField GaugeKmu(grid), Cmu(grid); GaugeLinkField GaugeKmu(grid), Cmu(grid);
StoutSmearing.BaseSmear(C, GaugeK); StoutSmearing.BaseSmear(C, GaugeK);
SigmaK = zero;
iLambda = zero;
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);
@ -81,14 +83,14 @@ GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
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,
@ -101,17 +103,17 @@ void set_iLambda(GaugeLinkField& iLambda,
GaugeLinkField unity(grid); GaugeLinkField unity(grid);
unity=1.0; unity=1.0;
LatticeReal u(grid), w(grid); LatticeComplex u(grid), w(grid);
LatticeComplex f0(grid), f1(grid), f2(grid); LatticeComplex f0(grid), f1(grid), f2(grid);
LatticeReal xi0(grid), xi1(grid), tmp(grid); LatticeComplex xi0(grid), xi1(grid), tmp(grid);
LatticeReal u2(grid), w2(grid), cosw(grid); LatticeComplex u2(grid), w2(grid), cosw(grid);
LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid); LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(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 LatticeUnitReal(grid); LatticeComplex LatticeUnitComplex(grid);
LatticeUnitReal = 1.0; LatticeUnitComplex = 1.0;
// Exponential // Exponential
iQ2 = iQ * iQ; iQ2 = iQ * iQ;
@ -121,44 +123,44 @@ void set_iLambda(GaugeLinkField& iLambda,
e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2; e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
// Getting B1, B2, Gamma and Lambda // Getting B1, B2, Gamma and Lambda
// simplify this part, reduntant calculations in set_fj
xi0 = StoutSmearing.func_xi0(w); xi0 = StoutSmearing.func_xi0(w);
xi1 = StoutSmearing.func_xi1(w); xi1 = StoutSmearing.func_xi1(w);
u2 = u * u; u2 = u * u;
w2 = w * w; w2 = w * w;
cosw = cos(w); cosw = cos(w);
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 = cos(u) - timesI(sin(u));
+ emiu * (toComplex(16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) + e2iu = cos(2.0*u) + timesI(sin(2.0*u));
timesI(toComplex(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0)));
r01 = (2.0*u + timesI(2.0*(u2-w2))) * e2iu
+ emiu * ((16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) +
timesI(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0));
r11 = (toComplex(2.0*LatticeUnitReal) + timesI(toComplex(4.0*u)))* e2iu r11 = (2.0*LatticeUnitComplex + timesI(4.0*u))* e2iu
+ emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) + + emiu * ((-2.0*cosw + (3.0*u2-w2)*xi0) +
timesI(toComplex(2.0*u*cosw + 6.0*u*xi0))); timesI((2.0*u*cosw + 6.0*u*xi0)));
r21 = 2.0*timesI(e2iu) r21 = 2.0*timesI(e2iu)
+ emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0))); + emiu * (-3.0*u*xi0 + timesI(cosw - 3.0*xi0));
r02 = -2.0 * e2iu + emiu * (toComplex(-8.0*u2*xi0) + r02 = -2.0 * e2iu + emiu * (-8.0*u2*xi0 +
timesI(toComplex(2.0*u*(cosw + xi0 + 3.0*u2*xi1)))); timesI(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 * (2.0*u*xi0 + timesI(-cosw - xi0 + 3.0*u2*xi1));
r22 = emiu * (toComplex(xi0) - timesI(toComplex(3.0*u*xi1))); r22 = emiu * (xi0 - timesI(3.0*u*xi1));
tmp = (2.0*(9.0*u2-w2)*(9.0*u2-w2)); fden = LatticeUnitComplex/(2.0*(9.0*u2-w2)*(9.0*u2-w2));
fden = toComplex(pow(tmp, -1.0)); // 1/tmp
b10 = 2.0 * u * r01 + (3.0* u2 - w2)*r02 - (30.0 * u2 + 2.0 * w2)*f0;
b11 = 2.0 * u * r11 + (3.0* u2 - w2)*r12 - (30.0 * u2 + 2.0 * w2)*f1;
b12 = 2.0 * u * r21 + (3.0* u2 - w2)*r22 - (30.0 * u2 + 2.0 * w2)*f2;
b10 = toComplex(2.0*u) * r01 + toComplex(3.0*u2 - w2)*r02 - toComplex(30.0*u2 + 2.0*w2)*f0; b20 = r01 - (3.0*u)*r02 - (24.0*u)*f0;
b11 = toComplex(2.0*u) * r11 + toComplex(3.0*u2 - w2)*r12 - toComplex(30.0*u2 + 2.0*w2)*f1; b21 = r11 - (3.0*u)*r12 - (24.0*u)*f1;
b12 = toComplex(2.0*u) * r21 + toComplex(3.0*u2 - w2)*r22 - toComplex(30.0*u2 + 2.0*w2)*f2; b22 = r21 - (3.0*u)*r22 - (24.0*u)*f2;
b20 = r01 - toComplex(3.0*u)*r02 - toComplex(24.0*u)*f0;
b21 = r11 - toComplex(3.0*u)*r12 - toComplex(24.0*u)*f1;
b22 = r21 - toComplex(3.0*u)*r22 - toComplex(24.0*u)*f2;
b10 *= fden; b10 *= fden;
b11 *= fden; b11 *= fden;
@ -166,7 +168,8 @@ void set_iLambda(GaugeLinkField& iLambda,
b20 *= fden; b20 *= fden;
b21 *= fden; b21 *= fden;
b22 *= fden; b22 *= fden;
B1 = b10*unity + timesMinusI(b11) * iQ - b12 * iQ2; B1 = b10*unity + timesMinusI(b11) * iQ - b12 * iQ2;
B2 = b20*unity + timesMinusI(b21) * iQ - b22 * iQ2; B2 = b20*unity + timesMinusI(b21) * iQ - b22 * iQ2;
USigmap = GaugeK * Sigmap; USigmap = GaugeK * Sigmap;
@ -180,58 +183,60 @@ void set_iLambda(GaugeLinkField& iLambda,
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(timesI(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
void set_GaugeField(GaugeField& U){ fill_smearedSet(U);} void set_GaugeField(GaugeField& U){ fill_smearedSet(U);}
//==================================================================== //====================================================================
void smeared_force(GaugeField& SigmaTilde) const{ void smeared_force(GaugeField& SigmaTilde) const{
if (smearingLevels > 0){
GaugeField force = SigmaTilde;//actually = U*SigmaTilde if (smearingLevels > 0){
GaugeLinkField tmp_mu(SigmaTilde._grid); GaugeField force(SigmaTilde._grid);
GaugeLinkField tmp_mu(SigmaTilde._grid);
for (int mu = 0; mu < Nd; mu++){ force = SigmaTilde;//actually = U*SigmaTilde
// to get SigmaTilde
tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu); for (int mu = 0; mu < Nd; mu++){
pokeLorentz(force, tmp_mu, mu); // to get just SigmaTilde
} tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
for(int ismr = smearingLevels - 1; ismr > 0; --ismr) pokeLorentz(force, tmp_mu, mu);
force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1)); }
force = AnalyticSmearedForce(force,*ThinLinks); for(int ismr = smearingLevels - 1; ismr > 0; --ismr)
force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1));
for (int mu = 0; mu < Nd; mu++){
tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); force = AnalyticSmearedForce(force,*ThinLinks);
pokeLorentz(SigmaTilde, tmp_mu, mu);
} 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 }// if smearingLevels = 0 do nothing
} }
//==================================================================== //====================================================================
@ -246,7 +251,7 @@ GaugeField& get_U(bool smeared=false) {
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 Usmr Plaq: " << impl_plaq<< std::endl;
return get_SmearedU(); return get_SmearedU();
} }

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

@ -12,7 +12,6 @@
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 Smear < Gimpl > * SmearBase; const Smear < Gimpl > * SmearBase;
public: public:
@ -22,58 +21,48 @@
static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3"); 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(){} //delete SmearBase...
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);
std::cout << "source matrix " << Umu << std::endl;
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);
//Debug check
GaugeLinkField check = adj(tmp) * tmp - 1.0;
std::cout << "check " << check << std::endl;
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,
const GaugeField& iLambda,
const GaugeField& Gauge) const{
SmearBase->derivative(SigmaTerm, iLambda, Gauge);
};
void derivative(GaugeField& SigmaTerm, void BaseSmear(GaugeField& C,
const GaugeField& iLambda, const GaugeField& U) const{
const GaugeField& Gauge) const{ SmearBase->smear(C, U);
SmearBase->derivative(SigmaTerm, iLambda, Gauge); };
};
void exponentiate_iQ(GaugeLinkField& e_iQ,
void BaseSmear(GaugeField& C, const GaugeLinkField& iQ) const{
const GaugeField& U) const{
SmearBase->smear(C, U);
};
void exponentiate_iQ(GaugeLinkField& e_iQ,
const GaugeLinkField& iQ) const{
// Put this outside // Put this outside
// only valid for SU(3) matrices // only valid for SU(3) matrices
@ -84,124 +73,48 @@
// the i sign is coming from outside // the i sign is coming from outside
// input matrix is anti-hermitian NOT hermitian // input matrix is anti-hermitian NOT hermitian
GridBase *grid = iQ._grid; GridBase *grid = iQ._grid;
GaugeLinkField unity(grid); GaugeLinkField unity(grid);
unity=1.0; unity=1.0;
GaugeLinkField iQ2(grid), iQ3(grid); GaugeLinkField iQ2(grid), iQ3(grid);
LatticeComplex u(grid), w(grid); LatticeComplex u(grid), w(grid);
LatticeComplex f0(grid), f1(grid), f2(grid); LatticeComplex f0(grid), f1(grid), f2(grid);
iQ2 = iQ * iQ; iQ2 = iQ * iQ;
iQ3 = iQ * iQ2; iQ3 = iQ * iQ2;
set_uw_complex(u, w, iQ2, iQ3); set_uw(u, w, iQ2, iQ3);
set_fj_complex(f0, f1, f2, u, w); set_fj(f0, f1, f2, u, w);
std::cout << "f0 " << f0 << std::endl; e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
std::cout << "f1 " << f1 << std::endl;
std::cout << "f2 " << f2 << std::endl;
std::cout << "iQ " << iQ << std::endl;
std::cout << "iQ2 " << iQ2 << std::endl;
e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
}; };
void set_uw(LatticeReal& u, LatticeReal& w, void set_uw(LatticeComplex& u, LatticeComplex& w,
GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{ GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
Real one_over_three = 1.0/3.0; Complex one_over_three = 1.0/3.0;
Real one_over_two = 1.0/2.0; Complex one_over_two = 1.0/2.0;
GridBase *grid = u._grid; GridBase *grid = u._grid;
LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid); LatticeComplex 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(real(timesMinusI(trace(iQ3)))) * one_over_three; //slow and 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);
theta = acos(c0/c0max);
u = sqrt(tmp) * cos( theta * one_over_three);
w = sqrt(c1) * sin( theta * one_over_three);
}
void set_uw_complex(LatticeComplex& u, LatticeComplex& w,
GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
Complex one_over_three = 1.0/3.0;
Complex one_over_two = 1.0/2.0;
GridBase *grid = u._grid;
LatticeComplex 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 = - real(timesMinusI(trace(iQ3))) * one_over_three; //temporary hack c0 = - real(timesMinusI(trace(iQ3))) * one_over_three; //temporary hack
c1 = - real(trace(iQ2)) * one_over_two; c1 = - real(trace(iQ2)) * one_over_two;
//Cayley Hamilton checks to machine precision, tested //Cayley Hamilton checks to machine precision, tested
std::cout << "c0 " << c0 << std::endl;
std::cout << "c1 " << c1 << std::endl;
tmp = c1 * one_over_three; tmp = c1 * one_over_three;
c0max = 2.0 * pow(tmp, 1.5); c0max = 2.0 * pow(tmp, 1.5);
std::cout << "c0max " << c0max << std::endl; theta = acos(c0/c0max)*one_over_three; // divide by three here, now leave as it is
LatticeComplex tempratio = c0/c0max; u = sqrt(tmp) * cos( theta );
std::cout << "ratio c0/c0max " << tempratio << std::endl; w = sqrt(c1) * sin( theta );
theta = acos(c0/c0max); // divide by three here, now leave as it is
std::cout << "theta " << theta << std::endl;
u = sqrt(tmp) * cos( theta * one_over_three);
w = sqrt(c1) * sin( theta * one_over_three);
std::cout << "u " << u << std::endl;
std::cout << "w " << w << std::endl;
} }
void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2, 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);
LatticeComplex fden(grid);
LatticeComplex h0(grid), h1(grid), h2(grid);
LatticeComplex e2iu(grid), emiu(grid), ixi0(grid), qt(grid);
xi0 = func_xi0(w);
u2 = u * u;
w2 = w * w;
cosw = cos(w);
ixi0 = timesI(toComplex(xi0));
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);
h1 = toComplex(2.0*u) * e2iu - emiu*( toComplex(2.0*u*cosw) -
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;
}
void set_fj_complex(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2,
const LatticeComplex& u, const LatticeComplex& w) const{ const LatticeComplex& u, const LatticeComplex& w) const{
GridBase *grid = u._grid; GridBase *grid = u._grid;
@ -212,47 +125,35 @@
LatticeComplex unity(grid); LatticeComplex unity(grid);
unity = 1.0; unity = 1.0;
xi0 = sin(w)/w;//func_xi0(w); xi0 = func_xi0(w);
std::cout << "xi0 " << xi0 << std::endl;
u2 = u * u; u2 = u * u;
std::cout << "u2 " << u2 << std::endl;
w2 = w * w; w2 = w * w;
std::cout << "w2 " << w2 << std::endl;
cosw = cos(w); cosw = cos(w);
std::cout << "cosw " << cosw << std::endl;
ixi0 = timesI(xi0); ixi0 = timesI(xi0);
emiu = cos(u) - timesI(sin(u)); emiu = cos(u) - timesI(sin(u));
e2iu = cos(2.0*u) + timesI(sin(2.0*u)); e2iu = cos(2.0*u) + timesI(sin(2.0*u));
std::cout << "emiu " << emiu << std::endl;
std::cout << "e2iu " << e2iu << std::endl;
h0 = e2iu * (u2 - w2) + emiu * ( (8.0*u2*cosw) + (2.0*u*(3.0*u2 + w2)*ixi0)); h0 = e2iu * (u2 - w2) + emiu * ( (8.0*u2*cosw) + (2.0*u*(3.0*u2 + w2)*ixi0));
h1 = e2iu * (2.0 * u) - emiu * ( (2.0*u*cosw) - (3.0*u2-w2)*ixi0); h1 = e2iu * (2.0 * u) - emiu * ( (2.0*u*cosw) - (3.0*u2-w2)*ixi0);
h2 = e2iu - emiu * ( cosw + (3.0*u)*ixi0); h2 = e2iu - emiu * ( cosw + (3.0*u)*ixi0);
std::cout << "h0 " << h0 << std::endl;
std::cout << "h1 " << h1 << std::endl;
std::cout << "h2 " << h2 << std::endl;
fden = unity/(9.0*u2 - w2);// reals fden = unity/(9.0*u2 - w2);// reals
std::cout << "fden " << fden << std::endl;
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{ LatticeComplex func_xi0(const LatticeComplex& 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{ LatticeComplex func_xi1(const LatticeComplex& 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);