diff --git a/Grid/qcd/action/gauge/Photon.h b/Grid/qcd/action/gauge/Photon.h index 6ff697c5..3ce64919 100644 --- a/Grid/qcd/action/gauge/Photon.h +++ b/Grid/qcd/action/gauge/Photon.h @@ -30,8 +30,9 @@ namespace Grid{ namespace QCD{ + template - class QedGimpl + class QedGImpl { public: typedef S Simd; @@ -43,27 +44,27 @@ namespace QCD{ typedef iImplGaugeLink SiteLink; typedef iImplGaugeField SiteField; - typedef SiteField SiteComplex; + typedef SiteLink SiteComplex; typedef Lattice LinkField; typedef Lattice Field; typedef Field ComplexField; }; - typedef QedGimpl QedGimplR; + typedef QedGImpl QedGImplR; - template + template class Photon { public: - INHERIT_GIMPL_TYPES(Gimpl); + INHERIT_GIMPL_TYPES(GImpl); + typedef typename SiteGaugeLink::scalar_object ScalarSite; + typedef typename ScalarSite::scalar_type ScalarComplex; GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3); - GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2, qedInf, 3); + GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2); public: - Photon(Gauge gauge, ZmScheme zmScheme); - Photon(Gauge gauge, ZmScheme zmScheme, std::vector improvements); - Photon(Gauge gauge, ZmScheme zmScheme, Real G0); - Photon(Gauge gauge, ZmScheme zmScheme, std::vector improvements, Real G0); + Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme, std::vector improvement); + Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme); virtual ~Photon(void) = default; void FreePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); @@ -73,345 +74,252 @@ namespace QCD{ const GaugeLinkField &weight); void UnitField(GaugeField &out); private: - void infVolPropagator(GaugeLinkField &out); - void invKHatSquared(GaugeLinkField &out); + void makeSpatialNorm(LatticeInteger &spNrm); + void makeKHat(std::vector &khat); + void makeInvKHatSquared(GaugeLinkField &out); void zmSub(GaugeLinkField &out); + void transverseProjectSpatial(GaugeField &out); + void gaugeTransform(GaugeField &out); private: - Gauge gauge_; - ZmScheme zmScheme_; - std::vector improvement_; - Real G0_; + GridBase *grid_; + Gauge gauge_; + ZmScheme zmScheme_; + std::vector improvement_; }; - typedef Photon PhotonR; + typedef Photon PhotonR; - template - Photon::Photon(Gauge gauge, ZmScheme zmScheme) - : gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector()), - G0_(0.15493339023106021408483720810737508876916113364521) - {} - - template - Photon::Photon(Gauge gauge, ZmScheme zmScheme, + template + Photon::Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme, std::vector improvements) - : gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), - G0_(0.15493339023106021408483720810737508876916113364521) + : grid_(grid), gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements) {} - template - Photon::Photon(Gauge gauge, ZmScheme zmScheme, Real G0) - : gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector()), G0_(G0) + template + Photon::Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme) + : Photon(grid, gauge, zmScheme, std::vector()) {} - template - Photon::Photon(Gauge gauge, ZmScheme zmScheme, - std::vector improvements, Real G0) - : gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), G0_(G0) - {} - - template - void Photon::FreePropagator (const GaugeField &in,GaugeField &out) + template + void Photon::FreePropagator(const GaugeField &in, GaugeField &out) { - FFT theFFT(in._grid); + FFT theFFT(dynamic_cast(grid_)); + GaugeField in_k(grid_); + GaugeField prop_k(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); + 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::infVolPropagator(GaugeLinkField &out) + template + void Photon::makeSpatialNorm(LatticeInteger &spNrm) { - auto *grid = dynamic_cast(out._grid); - LatticeReal xmu(grid); - GaugeLinkField one(grid); - const unsigned int nd = grid->_ndimension; - std::vector &l = grid->_fdimensions; - std::vector x0(nd,0); - TComplex Tone = Complex(1.0,0.0); - TComplex Tzero = Complex(G0_,0.0); - FFT fft(grid); - - one = Complex(1.0,0.0); - out = zero; - for(int mu = 0; mu < nd; mu++) + LatticeInteger coor(grid_); + std::vector l = grid_->FullDimensions(); + + spNrm = zero; + for(int mu = 0; mu < grid_->Nd() - 1; mu++) { - LatticeCoordinate(xmu,mu); - Real lo2 = l[mu]/2.0; - xmu = where(xmu < lo2, xmu, xmu-double(l[mu])); - out = out + toComplex(4*M_PI*M_PI*xmu*xmu); + LatticeCoordinate(coor, mu); + coor = where(coor < Integer(l[mu]/2), coor, coor - Integer(l[mu])); + spNrm = spNrm + coor*coor; } - pokeSite(Tone, out, x0); - out = one/out; - pokeSite(Tzero, out, x0); - fft.FFT_all_dim(out, out, FFT::forward); } - - template - void Photon::invKHatSquared(GaugeLinkField &out) + + template + void Photon::makeKHat(std::vector &khat) { - GridBase *grid = out._grid; - GaugeLinkField kmu(grid), one(grid); - const unsigned int nd = grid->_ndimension; - std::vector &l = grid->_fdimensions; - std::vector zm(nd,0); - TComplex Tone = Complex(1.0,0.0); - TComplex Tzero= Complex(0.0,0.0); - - one = Complex(1.0,0.0); - out = zero; - for(int mu = 0; mu < nd; mu++) + const unsigned int nd = grid_->Nd(); + std::vector l = grid_->FullDimensions(); + Complex ci(0., 1.); + + khat.resize(nd, grid_); + for (unsigned int mu = 0; mu < nd; ++mu) { Real twoPiL = M_PI*2./l[mu]; - - LatticeCoordinate(kmu,mu); - kmu = 2.*sin(.5*twoPiL*kmu); - out = out + kmu*kmu; + + LatticeCoordinate(khat[mu], mu); + khat[mu] = exp(0.5*ci*khat[mu])*2.*sin(.5*twoPiL*khat[mu]); } - pokeSite(Tone, out, zm); - out = one/out; - pokeSite(Tzero, out, zm); + } + + template + void Photon::makeInvKHatSquared(GaugeLinkField &out) + { + std::vector khat; + GaugeLinkField lone(grid_); + const unsigned int nd = grid_->Nd(); + std::vector zm(nd, 0); + ScalarSite one = ScalarComplex(1., 0.), z = ScalarComplex(0., 0.); + + out = zero; + makeKHat(khat); + for(int mu = 0; mu < nd; mu++) + { + out = out + khat[mu]*conjugate(khat[mu]); + } + lone = ScalarComplex(1., 0.); + pokeSite(one, out, zm); + out = lone/out; + pokeSite(z, out, zm); } - template - void Photon::zmSub(GaugeLinkField &out) + template + void Photon::zmSub(GaugeLinkField &out) { - GridBase *grid = out._grid; - const unsigned int nd = grid->_ndimension; - std::vector &l = grid->_fdimensions; - switch (zmScheme_) { case ZmScheme::qedTL: { - std::vector zm(nd,0); - TComplex Tzero = Complex(0.0,0.0); - - pokeSite(Tzero, out, zm); + std::vector zm(grid_->Nd(), 0); + ScalarSite z = ScalarComplex(0., 0.); + pokeSite(z, 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); - coor = where(coor < Integer(l[d]/2), coor, coor-Integer(l[d])); - spNrm = spNrm + coor*coor; - } - out = where(spNrm == Integer(0), 0.*out, out); + LatticeInteger spNrm(grid_); - // IR improvement + makeSpatialNorm(spNrm); + out = where(spNrm == Integer(0), 0.*out, out); for(int i = 0; i < improvement_.size(); i++) { - Real f = sqrt(improvement_[i]+1); - out = where(spNrm == Integer(i+1), f*out, out); + Real f = sqrt(improvement_[i] + 1); + out = where(spNrm == Integer(i + 1), f*out, out); } } default: + assert(0); break; } } - template - void Photon::MomentumSpacePropagator(const GaugeField &in, - GaugeField &out) + template + void Photon::transverseProjectSpatial(GaugeField &out) { - GridBase *grid = out._grid; - LatticeComplex momProp(grid); - - switch (zmScheme_) + const unsigned int nd = grid_->Nd(); + GaugeLinkField invKHat(grid_), cst(grid_); + LatticeInteger spNrm(grid_); + std::vector khat, a(nd, grid_), aProj(nd, grid_); + + invKHat = zero; + makeSpatialNorm(spNrm); + makeKHat(khat); + for (unsigned int mu = 0; mu < nd; ++mu) { - case ZmScheme::qedTL: - case ZmScheme::qedL: + a[mu] = peekLorentz(out, mu); + if (mu < nd - 1) { - invKHatSquared(momProp); - zmSub(momProp); - break; + invKHat += khat[mu]*conjugate(khat[mu]); } - case ZmScheme::qedInf: + } + cst = ScalarComplex(1., 0.); + invKHat = where(spNrm == Integer(0), cst, invKHat); + invKHat = cst/invKHat; + cst = zero; + invKHat = where(spNrm == Integer(0), cst, invKHat); + for (unsigned int mu = 0; mu < nd; ++mu) + { + aProj[mu] = a[mu]; + for (unsigned int nu = 0; nu < nd - 1; ++nu) { - infVolPropagator(momProp); - break; + aProj[mu] -= invKHat*khat[mu]*conjugate(khat[nu])*a[nu]; } + pokeLorentz(out, aProj[mu], mu); + } + } + + template + void Photon::gaugeTransform(GaugeField &out) + { + switch (gauge_) + { + case Gauge::feynman: + break; + case Gauge::coulomb: + transverseProjectSpatial(out); + break; + case Gauge::landau: + assert(0); + break; default: + assert(0); break; } + } + + template + void Photon::MomentumSpacePropagator(const GaugeField &in, + GaugeField &out) + { + LatticeComplex momProp(grid_); + + makeInvKHatSquared(momProp); + zmSub(momProp); out = in*momProp; } - template - void Photon::StochasticWeight(GaugeLinkField &weight) + template + void Photon::StochasticWeight(GaugeLinkField &weight) { - auto *grid = dynamic_cast(weight._grid); - const unsigned int nd = grid->_ndimension; - std::vector latt_size = grid->_fdimensions; - - switch (zmScheme_) + const unsigned int nd = grid_->Nd(); + std::vector l = grid_->FullDimensions(); + Integer vol = 1; + + for(unsigned int mu = 0; mu < nd; mu++) { - case ZmScheme::qedTL: - case ZmScheme::qedL: - { - Integer vol = 1; - for(int d = 0; d < nd; d++) - { - vol = vol * latt_size[d]; - } - invKHatSquared(weight); - weight = sqrt(vol)*sqrt(weight); - zmSub(weight); - break; - } - case ZmScheme::qedInf: - { - infVolPropagator(weight); - weight = sqrt(real(weight)); - break; - } - default: - break; + vol = vol*l[mu]; } + makeInvKHatSquared(weight); + weight = sqrt(vol)*sqrt(weight); + zmSub(weight); } - template - void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng) { - auto *grid = dynamic_cast(out._grid); - GaugeLinkField weight(grid); + GaugeLinkField weight(grid_); StochasticWeight(weight); StochasticField(out, rng, weight); } - template - void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng, + template + void Photon::StochasticField(GaugeField &out, GridParallelRNG &rng, const GaugeLinkField &weight) { - auto *grid = dynamic_cast(out._grid); - const unsigned int nd = grid->_ndimension; - GaugeLinkField r(grid); - GaugeField aTilde(grid); - FFT fft(grid); + const unsigned int nd = grid_->Nd(); + GaugeLinkField r(grid_); + GaugeField aTilde(grid_); + FFT fft(dynamic_cast(grid_)); - switch (zmScheme_) + for(unsigned int mu = 0; mu < nd; mu++) { - case ZmScheme::qedTL: - case ZmScheme::qedL: - { - for(int mu = 0; mu < nd; mu++) - { - gaussian(rng, r); - r = weight*r; - 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; + gaussian(rng, r); + r = weight*r; + pokeLorentz(aTilde, r, mu); } - + gaugeTransform(aTilde); fft.FFT_all_dim(out, aTilde, FFT::backward); - out = real(out); } - template - void Photon::UnitField(GaugeField &out) + template + void Photon::UnitField(GaugeField &out) { - auto *grid = dynamic_cast(out._grid); - const unsigned int nd = grid->_ndimension; - GaugeLinkField r(grid); + const unsigned int nd = grid_->Nd(); + GaugeLinkField r(grid_); - r = Complex(1.0,0.0); - - for(int mu = 0; mu < nd; mu++) + r = ScalarComplex(1., 0.); + for(unsigned int mu = 0; mu < nd; mu++) { pokeLorentz(out, r, mu); } - out = real(out); } -// 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 diff --git a/Hadrons/Modules/MGauge/StochEm.cc b/Hadrons/Modules/MGauge/StochEm.cc index 6f8bf55e..574387e4 100644 --- a/Hadrons/Modules/MGauge/StochEm.cc +++ b/Hadrons/Modules/MGauge/StochEm.cc @@ -70,7 +70,7 @@ void TStochEm::execute(void) LOG(Message) << "Generating stochastic EM potential..." << std::endl; std::vector improvements = strToVec(par().improvement); - PhotonR photon(par().gauge, par().zmScheme, improvements, par().G0_qedInf); + PhotonR photon(envGetGrid(EmField), par().gauge, par().zmScheme, improvements); auto &a = envGet(EmField, getName()); auto &w = envGet(EmComp, "_" + getName() + "_weight"); diff --git a/Hadrons/Modules/MGauge/StochEm.hpp b/Hadrons/Modules/MGauge/StochEm.hpp index a3f8cc96..b549387b 100644 --- a/Hadrons/Modules/MGauge/StochEm.hpp +++ b/Hadrons/Modules/MGauge/StochEm.hpp @@ -47,8 +47,7 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar, PhotonR::Gauge, gauge, PhotonR::ZmScheme, zmScheme, - std::string, improvement, - Real, G0_qedInf); + std::string, improvement); }; class TStochEm: public Module diff --git a/Hadrons/Modules/MGauge/UnitEm.cc b/Hadrons/Modules/MGauge/UnitEm.cc index d2ecad5e..97da8224 100644 --- a/Hadrons/Modules/MGauge/UnitEm.cc +++ b/Hadrons/Modules/MGauge/UnitEm.cc @@ -62,7 +62,7 @@ void TUnitEm::setup(void) // execution /////////////////////////////////////////////////////////////////// void TUnitEm::execute(void) { - PhotonR photon(0, 0); // Just chose arbitrary input values here + PhotonR photon(envGetGrid(EmField), 0, 0); // Just chose arbitrary input values here auto &a = envGet(EmField, getName()); LOG(Message) << "Generating unit EM potential..." << std::endl; photon.UnitField(a);