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

first (dirty) implementation of Feynman stoctachtic EM field

Browse Source
This commit is contained in:
Antonin Portelli 2016-10-21 13:10:13 +01:00
parent 997fd882ff
commit 63d219498b
1 changed files with 226 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,126 +1,234 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/gauge/Photon.h Source file: ./lib/qcd/action/gauge/Photon.h
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef QCD_PHOTON_ACTION_H #ifndef QCD_PHOTON_ACTION_H
#define QCD_PHOTON_ACTION_H #define QCD_PHOTON_ACTION_H
namespace Grid{ namespace Grid{
namespace QCD{ namespace QCD{
template<class Gimpl> template<class Gimpl>
class Photon { class Photon
{
public:
INHERIT_GIMPL_TYPES(Gimpl);
enum class Gauge {Feynman, Coulomb, Landau};
enum class ZmScheme {QedL, QedTL};
public:
Photon(Gauge gauge, ZmScheme zmScheme);
virtual ~Photon(void) = default;
void FreePropagator(const GaugeField &in, GaugeField &out);
void MomentumSpacePropagator(const GaugeField &in, GaugeField &out);
void StochasticField(GaugeField &out, GridParallelRNG &rng);
private:
void kHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out);
private:
Gauge gauge_;
ZmScheme zmScheme_;
};
public: template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme)
: gauge_(gauge), zmScheme_(zmScheme)
{}
INHERIT_GIMPL_TYPES(Gimpl); template<class Gimpl>
void Photon<Gimpl>::FreePropagator (const GaugeField &in,GaugeField &out)
{
FFT theFFT(in._grid);
enum PhotonType { FEYNMAN_L, FEYNMAN_TL }; GaugeField in_k(in._grid);
GaugeField prop_k(in._grid);
PhotonType GaugeBC; theFFT.FFT_all_dim(in_k,in,FFT::forward);
MomentumSpacePropagator(prop_k,in_k);
theFFT.FFT_all_dim(out,prop_k,FFT::backward);
}
Photon(PhotonType _GaugeBC) : GaugeBC(_GaugeBC){}; template<class Gimpl>
void Photon<Gimpl>::kHatSquared(GaugeLinkField &out)
{
GridBase *grid = out._grid;
GaugeLinkField kmu(grid);
const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
void FreePropagator (const GaugeField &in,GaugeField &out) { out = zero;
FFT theFFT((GridCartesian *) in._grid); for(int mu = 0; mu < nd; mu++)
{
RealD twoPiL = M_PI*2./l[mu];
GaugeField in_k(in._grid); LatticeCoordinate(kmu,mu);
GaugeField prop_k(in._grid); kmu = 2.*sin(.5*twoPiL*kmu);
out = out + kmu*kmu;
}
}
theFFT.FFT_all_dim(in_k,in,FFT::forward); template<class Gimpl>
MomentumSpacePropagator(prop_k,in_k); void Photon<Gimpl>::zmSub(GaugeLinkField &out)
theFFT.FFT_all_dim(out,prop_k,FFT::backward); {
GridBase *grid = out._grid;
const unsigned int nd = grid->_ndimension;
switch (zmScheme_)
{
case ZmScheme::QedTL:
{
std::vector<int> zm(nd,0);
TComplexD Tzero = ComplexD(0.0,0.0);
pokeSite(Tzero, out, zm);
break;
} }
void MomentumSpacePropagator(GaugeField &out,const GaugeField &in) { case ZmScheme::QedL:
if ( GaugeBC == FEYNMAN_L ) { {
FeynmanGaugeMomentumSpacePropagator_L(out,in); LatticeInteger spNrm(grid), coor(grid);
} else if ( GaugeBC == FEYNMAN_TL ) { GaugeLinkField z(grid);
FeynmanGaugeMomentumSpacePropagator_TL(out,in);
} else { // No coulomb yet spNrm = zero;
assert(0); for(int d = 0; d < grid->_ndimension - 1; d++)
} {
LatticeCoordinate(coor,d);
spNrm = spNrm + coor*coor;
}
out = where(spNrm == 0, 0.*out, out);
break;
} }
void FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, const GaugeField &in) { default:
break;
}
}
FeynmanGaugeMomentumSpacePropagator_TL(out,in); template<class Gimpl>
void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in,
GaugeField &out)
{
GridBase *grid = out._grid;
const unsigned int nd = grid->_ndimension;
LatticeComplex k2Inv(grid), one(grid);
GridBase *grid = out._grid; one = ComplexD(1.0,0.0);
LatticeInteger coor(grid); kHatSquared(k2Inv);
GaugeField zz(grid); zz=zero;
// xyzt std::vector<int> zm(nd,0);
for(int d = 0; d < grid->_ndimension-1;d++){ TComplexD Tone = ComplexD(1.0,0.0);
LatticeCoordinate(coor,d); TComplexD Tzero= ComplexD(0.0,0.0);
out = where(coor==Integer(0),zz,out);
}
}
void FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out, const GaugeField &in) { pokeSite(Tone, k2Inv, zm);
k2Inv = one/k2Inv;
pokeSite(Tzero, k2Inv,zm);
// what type LatticeComplex out = in*k2Inv;
GridBase *grid = out._grid; }
int nd = grid->_ndimension;
typedef typename GaugeField::vector_type vector_type; template<class Gimpl>
typedef typename GaugeField::scalar_type ScalComplex; void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng)
typedef Lattice<iSinglet<vector_type> > LatComplex; {
auto *grid = out._grid;
const unsigned int nd = grid->_ndimension;
GaugeLinkField k2(grid), r(grid);
GaugeField aTilde(grid);
FFT fft(grid);
std::vector<int> latt_size = grid->_fdimensions; kHatSquared(k2);
zmSub(k2);
LatComplex denom(grid); denom= zero; for(int mu = 0; mu < nd; mu++)
LatComplex one(grid); one = ScalComplex(1.0,0.0); {
LatComplex kmu(grid); gaussian(rng, r);
r = k2*r;
ScalComplex ci(0.0,1.0); pokeLorentz(aTilde, r, mu);
// momphase = n * 2pi / L }
for(int mu=0;mu<Nd;mu++) { fft.FFT_all_dim(out, aTilde, FFT::backward);
}
LatticeCoordinate(kmu,mu); // template<class Gimpl>
// void Photon<Gimpl>::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out,
RealD TwoPiL = M_PI * 2.0/ latt_size[mu]; // const GaugeField &in)
// {
kmu = TwoPiL * kmu ; //
// FeynmanGaugeMomentumSpacePropagator_TL(out,in);
denom = denom + 4.0*sin(kmu*0.5)*sin(kmu*0.5); // Wilson term //
} // GridBase *grid = out._grid;
std::vector<int> zero_mode(nd,0); // LatticeInteger coor(grid);
TComplexD Tone = ComplexD(1.0,0.0); // GaugeField zz(grid); zz=zero;
TComplexD Tzero= ComplexD(0.0,0.0); //
// // xyzt
pokeSite(Tone,denom,zero_mode); // for(int d = 0; d < grid->_ndimension-1;d++){
// LatticeCoordinate(coor,d);
denom= one/denom; // out = where(coor==Integer(0),zz,out);
// }
pokeSite(Tzero,denom,zero_mode); // }
//
out = zero; // template<class Gimpl>
out = in*denom; // void Photon<Gimpl>::FeynmanGaugeMomentumSpacePropagator_TL(GaugeField &out,
}; // const GaugeField &in)
// {
}; //
// // what type LatticeComplex
// GridBase *grid = out._grid;
// int nd = grid->_ndimension;
//
// typedef typename GaugeField::vector_type vector_type;
// typedef typename GaugeField::scalar_type ScalComplex;
// typedef Lattice<iSinglet<vector_type> > LatComplex;
//
// std::vector<int> latt_size = grid->_fdimensions;
//
// LatComplex denom(grid); denom= zero;
// LatComplex one(grid); one = ScalComplex(1.0,0.0);
// LatComplex kmu(grid);
//
// ScalComplex ci(0.0,1.0);
// // momphase = n * 2pi / L
// for(int mu=0;mu<Nd;mu++) {
//
// LatticeCoordinate(kmu,mu);
//
// RealD TwoPiL = M_PI * 2.0/ latt_size[mu];
//
// kmu = TwoPiL * kmu ;
//
// denom = denom + 4.0*sin(kmu*0.5)*sin(kmu*0.5); // Wilson term
// }
// std::vector<int> zero_mode(nd,0);
// TComplexD Tone = ComplexD(1.0,0.0);
// TComplexD Tzero= ComplexD(0.0,0.0);
//
// pokeSite(Tone,denom,zero_mode);
//
// denom= one/denom;
//
// pokeSite(Tzero,denom,zero_mode);
//
// out = zero;
// out = in*denom;
// };
}} }}
#endif #endif