From 63d219498ba4a036d74c8520e7455c94543cff49 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 21 Oct 2016 13:10:13 +0100 Subject: [PATCH] first (dirty) implementation of Feynman stoctachtic EM field --- lib/qcd/action/gauge/Photon.h | 344 ++++++++++++++++++++++------------ 1 file changed, 226 insertions(+), 118 deletions(-) diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h index 29d40f8b..ea6f21b9 100644 --- a/lib/qcd/action/gauge/Photon.h +++ b/lib/qcd/action/gauge/Photon.h @@ -1,126 +1,234 @@ - /************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./lib/qcd/action/gauge/Photon.h - - Copyright (C) 2015 - -Author: Peter Boyle - - 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 - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - 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., - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - - See the full license in the file "LICENSE" in the top level distribution directory - *************************************************************************************/ - /* END LEGAL */ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/gauge/Photon.h + + Copyright (C) 2015 + + Author: Peter Boyle + + 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 + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + 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., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ #ifndef QCD_PHOTON_ACTION_H #define QCD_PHOTON_ACTION_H namespace Grid{ - namespace QCD{ +namespace QCD{ + + template + 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_; + }; - template - class Photon { - - public: - - INHERIT_GIMPL_TYPES(Gimpl); - - enum PhotonType { FEYNMAN_L, FEYNMAN_TL }; - - PhotonType GaugeBC; - - Photon(PhotonType _GaugeBC) : GaugeBC(_GaugeBC){}; - - void FreePropagator (const GaugeField &in,GaugeField &out) { - FFT theFFT((GridCartesian *) in._grid); - - GaugeField in_k(in._grid); - GaugeField prop_k(in._grid); - - theFFT.FFT_all_dim(in_k,in,FFT::forward); - MomentumSpacePropagator(prop_k,in_k); - theFFT.FFT_all_dim(out,prop_k,FFT::backward); - } - void MomentumSpacePropagator(GaugeField &out,const GaugeField &in) { - if ( GaugeBC == FEYNMAN_L ) { - FeynmanGaugeMomentumSpacePropagator_L(out,in); - } else if ( GaugeBC == FEYNMAN_TL ) { - FeynmanGaugeMomentumSpacePropagator_TL(out,in); - } else { // No coulomb yet - assert(0); - } - } - void FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, const GaugeField &in) { - - FeynmanGaugeMomentumSpacePropagator_TL(out,in); - - GridBase *grid = out._grid; - LatticeInteger coor(grid); - GaugeField zz(grid); zz=zero; - - // xyzt - for(int d = 0; d < grid->_ndimension-1;d++){ - LatticeCoordinate(coor,d); - out = where(coor==Integer(0),zz,out); - } - } - - void 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 > LatComplex; + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme) + : gauge_(gauge), zmScheme_(zmScheme) + {} + + template + void Photon::FreePropagator (const GaugeField &in,GaugeField &out) + { + FFT theFFT(in._grid); - std::vector 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 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; - }; - - }; - + GaugeField in_k(in._grid); + GaugeField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + MomentumSpacePropagator(prop_k,in_k); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + } + + template + void Photon::kHatSquared(GaugeLinkField &out) + { + GridBase *grid = out._grid; + GaugeLinkField kmu(grid); + const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; + + out = zero; + for(int mu = 0; mu < nd; mu++) + { + RealD twoPiL = M_PI*2./l[mu]; + + LatticeCoordinate(kmu,mu); + kmu = 2.*sin(.5*twoPiL*kmu); + out = out + kmu*kmu; + } + } + + template + void Photon::zmSub(GaugeLinkField &out) + { + GridBase *grid = out._grid; + const unsigned int nd = grid->_ndimension; + + switch (zmScheme_) + { + case ZmScheme::QedTL: + { + std::vector zm(nd,0); + TComplexD Tzero = ComplexD(0.0,0.0); + + pokeSite(Tzero, out, zm); + + break; + } + case ZmScheme::QedL: + { + LatticeInteger spNrm(grid), coor(grid); + GaugeLinkField z(grid); + + spNrm = zero; + for(int d = 0; d < grid->_ndimension - 1; d++) + { + LatticeCoordinate(coor,d); + spNrm = spNrm + coor*coor; + } + out = where(spNrm == 0, 0.*out, out); + + break; + } + default: + break; + } + } + + template + void Photon::MomentumSpacePropagator(const GaugeField &in, + GaugeField &out) + { + GridBase *grid = out._grid; + const unsigned int nd = grid->_ndimension; + LatticeComplex k2Inv(grid), one(grid); + + one = ComplexD(1.0,0.0); + kHatSquared(k2Inv); + + std::vector 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; + } + + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) + { + auto *grid = out._grid; + const unsigned int nd = grid->_ndimension; + GaugeLinkField k2(grid), r(grid); + GaugeField aTilde(grid); + FFT fft(grid); + + kHatSquared(k2); + zmSub(k2); + for(int mu = 0; mu < nd; mu++) + { + gaussian(rng, r); + r = k2*r; + pokeLorentz(aTilde, r, mu); + } + fft.FFT_all_dim(out, aTilde, FFT::backward); + } +// template +// void Photon::FeynmanGaugeMomentumSpacePropagator_L(GaugeField &out, +// const GaugeField &in) +// { +// +// FeynmanGaugeMomentumSpacePropagator_TL(out,in); +// +// GridBase *grid = out._grid; +// LatticeInteger coor(grid); +// GaugeField zz(grid); zz=zero; +// +// // xyzt +// for(int d = 0; d < grid->_ndimension-1;d++){ +// LatticeCoordinate(coor,d); +// out = where(coor==Integer(0),zz,out); +// } +// } +// +// template +// void Photon::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 > LatComplex; +// +// std::vector 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 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