portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 15:55:37 +00:00
Code Releases Activity

Added infinite-volume photon in Photon.h (not checked yet)

Browse Source
This commit is contained in:
James Harrison 2017-10-11 16:44:51 -04:00
parent b99622d9fb
commit a80e43dbcf
1 changed files with 96 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
lib/qcd/action/gauge/Photon.h
View File

@ -58,7 +58,7 @@ namespace QCD{
public: public:
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3); GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3);
GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2); GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2, qedInf, 3);
public: public:
Photon(Gauge gauge, ZmScheme zmScheme); Photon(Gauge gauge, ZmScheme zmScheme);
virtual ~Photon(void) = default; virtual ~Photon(void) = default;
@ -69,6 +69,7 @@ namespace QCD{
void StochasticField(GaugeField &out, GridParallelRNG &rng, void StochasticField(GaugeField &out, GridParallelRNG &rng,
const GaugeLinkField &weight); const GaugeLinkField &weight);
private: private:
void infVolPropagator(GaugeLinkField &out);
void invKHatSquared(GaugeLinkField &out); void invKHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out); void zmSub(GaugeLinkField &out);
private: private:
@ -96,6 +97,32 @@ namespace QCD{
theFFT.FFT_all_dim(out,prop_k,FFT::backward); theFFT.FFT_all_dim(out,prop_k,FFT::backward);
} }
template<class Gimpl>
void Photon<Gimpl>::infVolPropagator(GaugeLinkField &out)
{
GridBase *grid = out._grid;
GaugeLinkField xmu(grid), one(grid);
const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
std::vector<int> x0(nd,0);
TComplex Tone = Complex(1.0,0.0);
TComplex Tzero = Complex(0.0,0.0);
FFT fft(grid);
one = Complex(1.0,0.0);
out = zero;
for(int mu = 0; mu < nd; mu++)
{
LatticeCoordinate(xmu,mu);
xmu = where(xmu < latt_size[mu]/2, xmu, xmu-latt_size[mu]/2);
out = out + 4*M_PI*M_PI*xmu*xmu;
}
pokeSite(Tone, out, x0);
out = one/out;
pokeSite(Tzero, out, x0);
fft.FFT_all_dim(out, out, FFT::forward);
}
template<class Gimpl> template<class Gimpl>
void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out) void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out)
{ {
@ -164,12 +191,27 @@ namespace QCD{
GaugeField &out) GaugeField &out)
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
LatticeComplex k2Inv(grid); LatticeComplex momProp(grid);
invKHatSquared(k2Inv); switch (zmScheme_)
zmSub(k2Inv); {
case ZmScheme::qedTL:
case ZmScheme::qedL:
{
invKHatSquared(momProp);
zmSub(momProp);
break;
}
case ZmScheme::qedInf:
{
infVolPropagator(momProp);
break;
}
default:
break;
}
out = in*k2Inv; out = in*momProp;
} }
template<class Gimpl> template<class Gimpl>
@ -179,6 +221,11 @@ namespace QCD{
const unsigned int nd = grid->_ndimension; const unsigned int nd = grid->_ndimension;
std::vector<int> latt_size = grid->_fdimensions; std::vector<int> latt_size = grid->_fdimensions;
switch (zmScheme_)
{
case ZmScheme::qedTL:
case ZmScheme::qedL:
{
Integer vol = 1; Integer vol = 1;
for(int d = 0; d < nd; d++) for(int d = 0; d < nd; d++)
{ {
@ -187,6 +234,17 @@ namespace QCD{
invKHatSquared(weight); invKHatSquared(weight);
weight = sqrt(vol*real(weight)); weight = sqrt(vol*real(weight));
zmSub(weight); zmSub(weight);
break;
}
case ZmScheme::qedInf:
{
infVolPropagator(weight);
weight = sqrt(real(weight));
break;
}
default:
break;
}
} }
template<class Gimpl> template<class Gimpl>
@ -209,12 +267,34 @@ namespace QCD{
GaugeField aTilde(grid); GaugeField aTilde(grid);
FFT fft(grid); FFT fft(grid);
switch (zmScheme_)
{
case ZmScheme::qedTL:
case ZmScheme::qedL:
{
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
gaussian(rng, r); gaussian(rng, r);
r = weight*r; r = weight*r;
pokeLorentz(aTilde, r, mu); pokeLorentz(aTilde, r, mu);
} }
break;
}
case ZmScheme::qedInf:
{
Complex shift(1., 1.); // This needs to be a GaugeLink element?
for(int mu = 0; mu < nd; mu++)
{
bernoulli(rng, r);
r = weight*(2.*r - shift);
pokeLorentz(aTilde, r, mu);
}
break;
}
default:
break;
}
fft.FFT_all_dim(out, aTilde, FFT::backward); fft.FFT_all_dim(out, aTilde, FFT::backward);
out = real(out); out = real(out);