portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Photon code fix

Browse Source
This commit is contained in:
Antonin Portelli 2018-12-14 17:36:38 +00:00
parent 856476a890
commit c4b472176c
1 changed files with 11 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Grid/qcd/action/gauge/Photon.h
View File

@ -137,10 +137,10 @@ namespace QCD{
khat.resize(nd, grid_); khat.resize(nd, grid_);
for (unsigned int mu = 0; mu < nd; ++mu) for (unsigned int mu = 0; mu < nd; ++mu)
{ {
Real twoPiL = M_PI*2./l[mu]; Real piL = M_PI/l[mu];
LatticeCoordinate(khat[mu], mu); LatticeCoordinate(khat[mu], mu);
khat[mu] = exp(0.5*ci*khat[mu])*2.*sin(.5*twoPiL*khat[mu]); khat[mu] = exp(piL*ci*khat[mu])*2.*sin(piL*khat[mu]);
} }
} }
@ -189,6 +189,7 @@ namespace QCD{
Real f = sqrt(improvement_[i] + 1); Real f = sqrt(improvement_[i] + 1);
out = where(spNrm == Integer(i + 1), f*out, out); out = where(spNrm == Integer(i + 1), f*out, out);
} }
break;
} }
default: default:
assert(0); assert(0);
@ -200,7 +201,7 @@ namespace QCD{
void Photon<GImpl>::transverseProjectSpatial(GaugeField &out) void Photon<GImpl>::transverseProjectSpatial(GaugeField &out)
{ {
const unsigned int nd = grid_->Nd(); const unsigned int nd = grid_->Nd();
GaugeLinkField invKHat(grid_), cst(grid_); GaugeLinkField invKHat(grid_), cst(grid_), spdiv(grid_);
LatticeInteger spNrm(grid_); LatticeInteger spNrm(grid_);
std::vector<GaugeLinkField> khat, a(nd, grid_), aProj(nd, grid_); std::vector<GaugeLinkField> khat, a(nd, grid_), aProj(nd, grid_);
@ -220,13 +221,15 @@ namespace QCD{
invKHat = cst/invKHat; invKHat = cst/invKHat;
cst = zero; cst = zero;
invKHat = where(spNrm == Integer(0), cst, invKHat); invKHat = where(spNrm == Integer(0), cst, invKHat);
spdiv = zero;
for (unsigned int nu = 0; nu < nd - 1; ++nu)
{
spdiv += conjugate(khat[nu])*a[nu];
}
spdiv *= invKHat;
for (unsigned int mu = 0; mu < nd; ++mu) for (unsigned int mu = 0; mu < nd; ++mu)
{ {
aProj[mu] = a[mu]; aProj[mu] = a[mu] - khat[mu]*spdiv;
for (unsigned int nu = 0; nu < nd - 1; ++nu)
{
aProj[mu] -= invKHat*khat[mu]*conjugate(khat[nu])*a[nu];
}
pokeLorentz(out, aProj[mu], mu); pokeLorentz(out, aProj[mu], mu);
} }
} }