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

QED: fix stochastic field

Browse Source
This commit is contained in:
Antonin Portelli 2016-10-21 14:41:08 +01:00
parent bd6a228af6
commit 462921e549
1 changed files with 18 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -45,7 +45,7 @@ namespace QCD{
void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out);
void StochasticField(GaugeField &out, GridParallelRNG &rng); void StochasticField(GaugeField &out, GridParallelRNG &rng);
private: private:
void kHatSquared(GaugeLinkField &out); void invKHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out); void zmSub(GaugeLinkField &out);
private: private:
Gauge gauge_; Gauge gauge_;
@ -71,13 +71,17 @@ namespace QCD{
} }
template<class Gimpl> template<class Gimpl>
void Photon<Gimpl>::kHatSquared(GaugeLinkField &out) void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out)
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
GaugeLinkField kmu(grid); GaugeLinkField kmu(grid), one(grid);
const unsigned int nd = grid->_ndimension; const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions; std::vector<int> &l = grid->_fdimensions;
std::vector<int> zm(nd,0);
TComplexD Tone = ComplexD(1.0,0.0);
TComplexD Tzero= ComplexD(0.0,0.0);
one = ComplexD(1.0,0.0);
out = zero; out = zero;
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
@ -87,6 +91,9 @@ namespace QCD{
kmu = 2.*sin(.5*twoPiL*kmu); kmu = 2.*sin(.5*twoPiL*kmu);
out = out + kmu*kmu; out = out + kmu*kmu;
} }
pokeSite(Tone, out, zm);
out = one/out;
pokeSite(Tzero, out,zm);
} }
template<class Gimpl> template<class Gimpl>
@ -131,19 +138,10 @@ namespace QCD{
GaugeField &out) GaugeField &out)
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
const unsigned int nd = grid->_ndimension; LatticeComplex k2Inv(grid);
LatticeComplex k2Inv(grid), one(grid);
one = ComplexD(1.0,0.0); invKHatSquared(k2Inv);
kHatSquared(k2Inv); zmSub(k2Inv);
std::vector<int> zm(nd,0);
TComplexD Tone = ComplexD(1.0,0.0);
TComplexD Tzero= ComplexD(0.0,0.0);
pokeSite(Tone, k2Inv, zm);
k2Inv = one/k2Inv;
pokeSite(Tzero, k2Inv,zm);
out = in*k2Inv; out = in*k2Inv;
} }
@ -153,16 +151,17 @@ namespace QCD{
{ {
auto *grid = out._grid; auto *grid = out._grid;
const unsigned int nd = grid->_ndimension; const unsigned int nd = grid->_ndimension;
GaugeLinkField k2(grid), r(grid); GaugeLinkField sqrtK2Inv(grid), r(grid);
GaugeField aTilde(grid); GaugeField aTilde(grid);
FFT fft(grid); FFT fft(grid);
kHatSquared(k2); invKHatSquared(sqrtK2Inv);
zmSub(k2); sqrtK2Inv = sqrt(real(sqrtK2Inv));
zmSub(sqrtK2Inv);
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
gaussian(rng, r); gaussian(rng, r);
r = k2*r; r = sqrtK2Inv*r;
pokeLorentz(aTilde, r, mu); pokeLorentz(aTilde, r, mu);
} }
fft.FFT_all_dim(out, aTilde, FFT::backward); fft.FFT_all_dim(out, aTilde, FFT::backward);