portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-13 01:05:36 +00:00
Code Releases Activity

stochastic QED: function to cache 1/sqrt(khat^2)

Browse Source
This commit is contained in:
Antonin Portelli 2016-12-22 00:29:19 +01:00
parent 149a46b92c
commit 17b3a10d46
1 changed files with 34 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -44,7 +44,10 @@ namespace QCD{
virtual ~Photon(void) = default; virtual ~Photon(void) = default;
void FreePropagator(const GaugeField &in, GaugeField &out); void FreePropagator(const GaugeField &in, GaugeField &out);
void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out);
void StochasticWeight(GaugeLinkField &weight);
void StochasticField(GaugeField &out, GridParallelRNG &rng); void StochasticField(GaugeField &out, GridParallelRNG &rng);
void StochasticField(GaugeField &out, GridParallelRNG &rng,
const GaugeLinkField &weight);
private: private:
void invKHatSquared(GaugeLinkField &out); void invKHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out); void zmSub(GaugeLinkField &out);
@ -148,32 +151,50 @@ namespace QCD{
} }
template<class Gimpl> template<class Gimpl>
void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng) void Photon<Gimpl>::StochasticWeight(GaugeLinkField &weight)
{ {
auto *grid = dynamic_cast<GridCartesian *>(out._grid); auto *grid = dynamic_cast<GridCartesian *>(weight._grid);
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;
GaugeLinkField sqrtK2Inv(grid), r(grid);
GaugeField aTilde(grid);
FFT fft(grid);
Integer vol = 1; Integer vol = 1;
for(int d = 0; d < nd; d++) for(int d = 0; d < nd; d++)
{ {
vol = vol * latt_size[d]; vol = vol * latt_size[d];
} }
invKHatSquared(weight);
invKHatSquared(sqrtK2Inv); weight = sqrt(vol*real(weight));
sqrtK2Inv = sqrt(vol*real(sqrtK2Inv)); zmSub(weight);
zmSub(sqrtK2Inv); }
template<class Gimpl>
void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng)
{
auto *grid = dynamic_cast<GridCartesian *>(out._grid);
GaugeLinkField weight(grid);
StochasticWeight(weight);
StochasticField(out, rng, weight);
}
template<class Gimpl>
void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng,
const GaugeLinkField &weight)
{
auto *grid = dynamic_cast<GridCartesian *>(out._grid);
const unsigned int nd = grid->_ndimension;
GaugeLinkField r(grid);
GaugeField aTilde(grid);
FFT fft(grid);
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
gaussian(rng, r); gaussian(rng, r);
r = sqrtK2Inv*r; r = weight*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);
out = real(out); out = real(out);
} }
// template<class Gimpl> // template<class Gimpl>