portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-03 21:44:33 +00:00
Code Releases Activity

Smearing routines compile (still untested)

Browse Source
This commit is contained in:
Guido Cossu
2016-02-25 02:43:59 +09:00
parent a7251f28c7
commit 2d8bb356e3
8 changed files with 249 additions and 92 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/lattice/Lattice_ET.h
View File

@@ -213,6 +213,7 @@ inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3
GridUnopClass(UnaryToReal,toReal(a));
GridUnopClass(UnaryToComplex,toComplex(a));
GridUnopClass(UnaryTimesI,timesI(a));
GridUnopClass(UnaryTimesMinusI,timesMinusI(a));
GridUnopClass(UnaryAbs,abs(a));
GridUnopClass(UnarySqrt,sqrt(a));
GridUnopClass(UnaryRsqrt,rsqrt(a));
@@ -329,6 +330,7 @@ GRID_DEF_UNOP(imag,UnaryImag);
GRID_DEF_UNOP(toReal,UnaryToReal);
GRID_DEF_UNOP(toComplex,UnaryToComplex);
GRID_DEF_UNOP(timesI,UnaryTimesI);
GRID_DEF_UNOP(timesMinusI,UnaryTimesMinusI);
GRID_DEF_UNOP(abs ,UnaryAbs); //abs overloaded in cmath C++98; DON'T do the abs-fabs-dabs-labs thing
GRID_DEF_UNOP(sqrt ,UnarySqrt);
GRID_DEF_UNOP(rsqrt,UnaryRsqrt);

2
lib/qcd/action/gauge/GaugeImpl.h
View File

@@ -61,7 +61,7 @@ template<class Gimpl> class WilsonLoops;
typedef Lattice<SiteGaugeField> GaugeField;
// Move this elsewhere?
void AddGaugeLink(GaugeField& U, GaugeLinkField& W, int mu){ // U[mu] += W
static inline void AddGaugeLink(GaugeField& U, GaugeLinkField& W, int mu){ // U[mu] += W
PARALLEL_FOR_LOOP
for(auto ss=0;ss<U._grid->oSites();ss++){
U._odata[ss]._internal[mu] = U._odata[ss]._internal[mu] + W._odata[ss]._internal;

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

@@ -86,19 +86,19 @@ namespace Grid {
temp_Sigma = -rho_numu*staple*iLambda_nu;
//-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
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)
AddGaugeLink(SigmaTerm, temp_Sigma, mu);
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
sh_field = Cshift(iLambda_mu, nu, 1);
temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu);
//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
AddGaugeLink(SigmaTerm, temp_Sigma, mu);
Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
staple = zero;
sh_field = Cshift(U_nu, mu, 1);
@@ -110,7 +110,7 @@ namespace Grid {
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);
AddGaugeLink(SigmaTerm, sh_field, mu);
Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
}
}

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

@@ -24,34 +24,149 @@ namespace Grid {
template <class Gimpl>
class SmearedConfiguration {
public:
INHERIT_GIMPL_TYPES(Gimpl)
INHERIT_GIMPL_TYPES(Gimpl) ;
private:
const unsigned int smearingLevels;
Smear_Stout StoutSmearing;
Smear_Stout<Gimpl> StoutSmearing;
std::vector<GaugeField> SmearedSet;
// Member functions
void fill_smearedSet();
GaugeField AnalyticSmearedForce(const GaugeField&,
const GaugeField&) const;
const GaugeField& get_smeared_conf(int) const;
//====================================================================
void fill_smearedSet(){
GaugeField previous_u;
void set_iLambda(GaugeField& iLambda,
GaugeField& e_iQ,
const GaugeField& iQ,
const GaugeField& Sigmap,
const GaugeField& U)const;
std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
/* Check these types (do I need to pass iQ1,2 ? )
void set_uw(RealD& u, RealD& w,
const SUNmat& iQ1, const SUNmat& iQ2)const ;
void set_fj(ComplexD& f0, ComplexD& f1,
CompledD& f2, const RealD& u,
const RealD& w)const;
*/
previous_u = *ThinLinks;
for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
previous_u = SmearedSet[smearLvl];
}
RealD func_xi0(RealD w)const;
RealD func_xi1(RealD w)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);
GaugeLinkField iQ(grid), e_iQ(grid);
GaugeLinkField SigmaKPrime_mu(grid);
GaugeLinkField GaugeKmu(grid), Cmu(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);
}
StoutSmearing.derivative(SigmaK, iLambda, GaugeK);
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,
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);
unity=1.0;
LatticeReal u(grid), w(grid);
LatticeComplex f0(grid), f1(grid), f2(grid);
LatticeReal xi0(grid), xi1(grid), fden(grid), tmp(grid);
LatticeReal u2(grid), w2(grid), cosw(grid);
LatticeComplex emiu(grid), e2iu(grid), qt(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);
unitReal = 1.0;
// Exponential
iQ2 = iQ * iQ;
iQ3 = iQ * iQ2;
StoutSmearing.set_uw(u,w,iQ2,iQ3);
StoutSmearing.set_fj(f0,f1,f2,u,w);
e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
// Getting B1, B2, Gamma and Lambda
xi0 = StoutSmearing.func_xi0(w);
xi1 = StoutSmearing.func_xi1(w);
u2 = u * u;
w2 = w * w;
cosw = cos(w);
emiu = toComplex(cos(u)) - timesI(toComplex(u));
e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(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)));
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)));
r21 = timesI(toComplex(2.0*unitReal)) * 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))));
r12 = emiu * (toComplex(2.0*u*xi0) + timesI(toComplex(-cosw - xi0 + 3.0*u2*xi1)));
r22 = emiu * (toComplex(xi0) - timesI(toComplex(3.0*u*xi1)));
tmp = (2.0*(9.0*u2-w2)*(9.0*u2-w2));
fden = toComplex(pow(tmp, -1.0)); // 1/tmp
b10 = toComplex(2.0*u) * r01 + toComplex(3.0*u2 - w2)*r02 - toComplex(30.0*u2 + 2.0*w2)*f0;
b11 = toComplex(2.0*u) * r11 + toComplex(3.0*u2 - w2)*r12 - toComplex(30.0*u2 + 2.0*w2)*f1;
b12 = toComplex(2.0*u) * r21 + toComplex(3.0*u2 - w2)*r22 - toComplex(30.0*u2 + 2.0*w2)*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;
b11 *= fden;
b12 *= fden;
b20 *= fden;
b21 *= fden;
b22 *= fden;
B1 = b10*unity + timesMinusI(b11) * iQ - b12 * iQ2;
B2 = b20*unity + timesMinusI(b21) * iQ - b22 * iQ2;
USigmap = GaugeK * Sigmap;
tr1 = trace(USigmap*B1);
tr2 = trace(USigmap*B2);
GaugeLinkField QUS = timesMinusI(iQ) * USigmap;
GaugeLinkField USQ = USigmap * timesMinusI(iQ);
GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 +
f1 * USigmap + f2 * QUS + f2 * USQ;
iLambda = Ta(iGamma);
}
public:
GaugeField* ThinLinks; /*!< @brief Pointer to the thin
@@ -60,10 +175,10 @@ namespace Grid {
/*! @brief Standard constructor */
SmearedConfiguration(GridCartesian * UGrid,
unsigned int Nsmear,
Smear_Stout& Stout):
Smear_Stout<Gimpl>& Stout):
smearingLevels(Nsmear),
StoutSmearing(Stout),
ThinLinks(new GaugeField){
ThinLinks(new GaugeField(UGrid)){
for (unsigned int i=0; i< smearingLevels; ++i)
SmearedSet.push_back(*(new GaugeField(UGrid)));
}
@@ -76,7 +191,26 @@ namespace Grid {
ThinLinks(new GaugeField(UGrid)){}
void set_GaugeField(){ fill_smearedSet(); }
void smeared_force(GaugeField&) const;
void smeared_force(GaugeField& SigmaTilde) const{
GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid
GaugeLinkField tmp_mu(SigmaTilde._grid);
for (int mu = 0; mu < Nd; 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,*ThinLinks);
for (int mu = 0; mu < Nd; mu++){
tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
pokeLorentz(SigmaTilde, tmp_mu, mu);
}
}
GaugeField* get_SmearedU() const{
return const_cast<GaugeField*>(&(SmearedSet[smearingLevels-1]));
}

2
lib/qcd/smearing/Smearing.h
View File

@@ -4,6 +4,6 @@
#include <qcd/smearing/BaseSmearing.h>
#include <qcd/smearing/APEsmearing.h>
#include <qcd/smearing/StoutSmearing.h>
#include <qcd/smearing/GaugeConfiguration.h>
#endif

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

@@ -15,11 +15,7 @@ namespace Grid {
const std::vector<double> d_rho;
const Smear < Gimpl > * SmearBase;
LatticeReal func_xi0(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;
}
public:
INHERIT_GIMPL_TYPES(Gimpl)
@@ -32,20 +28,20 @@ namespace Grid {
~Smear_Stout(){}
void smear(GaugeField& u_smr,const GaugeField& U) const{
long double timing;
GaugeField u_tmp1, q_mu;
GaugeField C(U._grid);
GaugeLinkField tmp(U._grid), q_mu(U._grid), Umu(U._grid);
std::cout<< GridLogDebug << "Stout smearing started\n";
//Smear the configurations
SmearBase->smear(u_tmp1, U);
q_mu = Ta(u_tmp1*adj(u_tmp1)); // q_mu = Ta(Omega_mu)
exponentiate_iQ(u_tmp1, q_mu);
u_smr = u_tmp1*U;
SmearBase->smear(C, U);
for (int mu = 0; mu<Nd; mu++){
tmp = peekLorentz(C,mu);
Umu = peekLorentz(U,mu);
q_mu = Ta(tmp * adj(Umu)); // q_mu = Ta(Omega_mu)
exponentiate_iQ(tmp, q_mu);
pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu
}
std::cout<< GridLogDebug << "Stout smearing completed\n";
}
@@ -61,58 +57,74 @@ namespace Grid {
SmearBase->smear(C, U);
}
void exponentiate_iQ(GaugeField& e_iQ,
const GaugeField& iQ) const{
void exponentiate_iQ(GaugeLinkField& e_iQ,
const GaugeLinkField& iQ) const{
// Put this outside
// only valid for SU(3) matrices
// only one Lorentz direction at a time
GridBase *grid = iQ._grid;
Real one_over_three = 1.0/3.0;
Real one_over_two = 1.0/2.0;
GaugeLinkField unity(grid);
unity=1.0;
GaugeField unity;
GaugeLinkField Umu(iQ._grid);
Umu=1.0;
for(int mu=0;mu<Nd;mu++){
pokeLorentz(unity,Umu,mu);
}
GaugeField iQ2, iQ3;
LatticeReal c0(grid), c1(grid), c0max(grid), u_val(grid), tmp(grid);
LatticeReal w(grid), theta(grid), xi0(grid), u2(grid), w2(grid), cosw(grid);
LatticeComplex fden(grid);
LatticeComplex f0(grid), f1(grid), f2(grid), h0(grid), h1(grid), h2(grid);
LatticeComplex e2iu(grid), emiu(grid), ixi0(grid), qt(grid);
GaugeLinkField iQ2(grid), iQ3(grid);
LatticeReal u(grid), w(grid);
LatticeComplex f0(grid), f1(grid), f2(grid);
iQ2 = iQ * iQ;
iQ3 = iQ * iQ2;
c0 = - imag(trace(iQ3)) * one_over_three;
c1 = - real(trace(iQ2)) * one_over_two;
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{
Real one_over_three = 1.0/3.0;
Real one_over_two = 1.0/2.0;
GridBase *grid = u._grid;
LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
c0 = - toReal(imag(trace(iQ3))) * one_over_three;
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_val = sqrt(tmp) * cos( theta * one_over_three);
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{
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_val * u_val;
u2 = u * u;
w2 = w * w;
cosw = cos(w);
ixi0 = timesI(toComplex(xi0));
emiu = toComplex(cos(u_val)) - timesI(toComplex(u_val));
e2iu = toComplex(cos(2.0*u_val)) + timesI(toComplex(2.0*u_val));
emiu = toComplex(cos(u)) - timesI(toComplex(u));
e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u));
h0 = e2iu * toComplex(u2 - w2) + emiu *( toComplex(8.0*u2*cosw) +
toComplex(2.0*u_val*(3.0*u2 + w2))*ixi0);
toComplex(2.0*u*(3.0*u2 + w2))*ixi0);
h1 = toComplex(2.0*u_val) * e2iu - emiu*( toComplex(2.0*u_val*cosw) -
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_val)*ixi0);
h2 = e2iu - emiu * (toComplex(cosw) + toComplex(3.0*u)*ixi0);
tmp = 9.0*u2 - w2;
@@ -120,13 +132,20 @@ namespace Grid {
f0 = h0 * fden;
f1 = h1 * fden;
f2 = h2 * fden;
}
e_iQ = f0*unity + f1 * timesMinusI(iQ) - f2 * iQ2;
LatticeReal func_xi0(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(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);
}
};

4
lib/simd/Grid_vector_unops.h
View File

@@ -118,12 +118,12 @@ namespace Grid {
};
template<class scalar> struct RealFunctor {
scalar operator()(const scalar &a) const {
scalar operator()(const std::complex<scalar> &a) const {
return real(a);
}
};
template<class scalar> struct ImagFunctor {
scalar operator()(const scalar &a) const {
scalar operator()(const std::complex<scalar> &a) const {
return imag(a);
}
};

4
tests/Test_main.cc
View File

@@ -609,8 +609,10 @@ int main (int argc, char ** argv)
// Testing Smearing routine compilation
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation
Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing);
SmearedConfiguration< PeriodicGimplR > SmartConf(&Fine, 3, StoutSmearing);
std::cout<<GridLogMessage << sizeof(vComplexF) << std::endl;