portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

big cleanup of the Photon class + QED Coulomb gauge

Browse Source
This commit is contained in:
Antonin Portelli 2018-12-13 21:52:38 +00:00
parent c509bd3fe2
commit 856476a890
4 changed files with 176 additions and 269 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

438
Grid/qcd/action/gauge/Photon.h
View File

@ -30,8 +30,9 @@
namespace Grid{ namespace Grid{
namespace QCD{ namespace QCD{
template <class S> template <class S>
class QedGimpl class QedGImpl
{ {
public: public:
typedef S Simd; typedef S Simd;
@ -43,27 +44,27 @@ namespace QCD{
typedef iImplGaugeLink<Simd> SiteLink; typedef iImplGaugeLink<Simd> SiteLink;
typedef iImplGaugeField<Simd> SiteField; typedef iImplGaugeField<Simd> SiteField;
typedef SiteField SiteComplex; typedef SiteLink SiteComplex;
typedef Lattice<SiteLink> LinkField; typedef Lattice<SiteLink> LinkField;
typedef Lattice<SiteField> Field; typedef Lattice<SiteField> Field;
typedef Field ComplexField; typedef Field ComplexField;
}; };
typedef QedGimpl<vComplex> QedGimplR; typedef QedGImpl<vComplex> QedGImplR;
template<class Gimpl> template <class GImpl>
class Photon class Photon
{ {
public: 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(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: public:
Photon(Gauge gauge, ZmScheme zmScheme); Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme, std::vector<Real> improvement);
Photon(Gauge gauge, ZmScheme zmScheme, std::vector<Real> improvements); Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme);
Photon(Gauge gauge, ZmScheme zmScheme, Real G0);
Photon(Gauge gauge, ZmScheme zmScheme, std::vector<Real> improvements, Real G0);
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);
@ -73,345 +74,252 @@ namespace QCD{
const GaugeLinkField &weight); const GaugeLinkField &weight);
void UnitField(GaugeField &out); void UnitField(GaugeField &out);
private: private:
void infVolPropagator(GaugeLinkField &out); void makeSpatialNorm(LatticeInteger &spNrm);
void invKHatSquared(GaugeLinkField &out); void makeKHat(std::vector<GaugeLinkField> &khat);
void makeInvKHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out); void zmSub(GaugeLinkField &out);
void transverseProjectSpatial(GaugeField &out);
void gaugeTransform(GaugeField &out);
private: private:
Gauge gauge_; GridBase *grid_;
ZmScheme zmScheme_; Gauge gauge_;
std::vector<Real> improvement_; ZmScheme zmScheme_;
Real G0_; std::vector<Real> improvement_;
}; };
typedef Photon<QedGimplR> PhotonR; typedef Photon<QedGImplR> PhotonR;
template<class Gimpl> template<class GImpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme) Photon<GImpl>::Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme,
: gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector<Real>()),
G0_(0.15493339023106021408483720810737508876916113364521)
{}
template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme,
std::vector<Real> improvements) std::vector<Real> improvements)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), : grid_(grid), gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements)
G0_(0.15493339023106021408483720810737508876916113364521)
{} {}
template<class Gimpl> template<class GImpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme, Real G0) Photon<GImpl>::Photon(GridBase *grid, Gauge gauge, ZmScheme zmScheme)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector<Real>()), G0_(G0) : Photon(grid, gauge, zmScheme, std::vector<Real>())
{} {}
template<class Gimpl> template<class GImpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme, void Photon<GImpl>::FreePropagator(const GaugeField &in, GaugeField &out)
std::vector<Real> improvements, Real G0)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), G0_(G0)
{}
template<class Gimpl>
void Photon<Gimpl>::FreePropagator (const GaugeField &in,GaugeField &out)
{ {
FFT theFFT(in._grid); FFT theFFT(dynamic_cast<GridCartesian *>(grid_));
GaugeField in_k(grid_);
GaugeField prop_k(grid_);
GaugeField in_k(in._grid); theFFT.FFT_all_dim(in_k, in, FFT::forward);
GaugeField prop_k(in._grid); 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<class Gimpl> template<class GImpl>
void Photon<Gimpl>::infVolPropagator(GaugeLinkField &out) void Photon<GImpl>::makeSpatialNorm(LatticeInteger &spNrm)
{ {
auto *grid = dynamic_cast<GridCartesian *>(out._grid); LatticeInteger coor(grid_);
LatticeReal xmu(grid); std::vector<int> l = grid_->FullDimensions();
GaugeLinkField one(grid);
const unsigned int nd = grid->_ndimension; spNrm = zero;
std::vector<int> &l = grid->_fdimensions; for(int mu = 0; mu < grid_->Nd() - 1; mu++)
std::vector<int> 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++)
{ {
LatticeCoordinate(xmu,mu); LatticeCoordinate(coor, mu);
Real lo2 = l[mu]/2.0; coor = where(coor < Integer(l[mu]/2), coor, coor - Integer(l[mu]));
xmu = where(xmu < lo2, xmu, xmu-double(l[mu])); spNrm = spNrm + coor*coor;
out = out + toComplex(4*M_PI*M_PI*xmu*xmu);
} }
pokeSite(Tone, out, x0);
out = one/out;
pokeSite(Tzero, out, x0);
fft.FFT_all_dim(out, out, FFT::forward);
} }
template<class Gimpl> template<class GImpl>
void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out) void Photon<GImpl>::makeKHat(std::vector<GaugeLinkField> &khat)
{ {
GridBase *grid = out._grid; const unsigned int nd = grid_->Nd();
GaugeLinkField kmu(grid), one(grid); std::vector<int> l = grid_->FullDimensions();
const unsigned int nd = grid->_ndimension; Complex ci(0., 1.);
std::vector<int> &l = grid->_fdimensions;
std::vector<int> zm(nd,0); khat.resize(nd, grid_);
TComplex Tone = Complex(1.0,0.0); for (unsigned int mu = 0; mu < nd; ++mu)
TComplex Tzero= Complex(0.0,0.0);
one = Complex(1.0,0.0);
out = zero;
for(int mu = 0; mu < nd; mu++)
{ {
Real twoPiL = M_PI*2./l[mu]; Real twoPiL = M_PI*2./l[mu];
LatticeCoordinate(kmu,mu); LatticeCoordinate(khat[mu], mu);
kmu = 2.*sin(.5*twoPiL*kmu); khat[mu] = exp(0.5*ci*khat[mu])*2.*sin(.5*twoPiL*khat[mu]);
out = out + kmu*kmu;
} }
pokeSite(Tone, out, zm); }
out = one/out;
pokeSite(Tzero, out, zm); template<class GImpl>
void Photon<GImpl>::makeInvKHatSquared(GaugeLinkField &out)
{
std::vector<GaugeLinkField> khat;
GaugeLinkField lone(grid_);
const unsigned int nd = grid_->Nd();
std::vector<int> 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<class Gimpl> template<class GImpl>
void Photon<Gimpl>::zmSub(GaugeLinkField &out) void Photon<GImpl>::zmSub(GaugeLinkField &out)
{ {
GridBase *grid = out._grid;
const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
switch (zmScheme_) switch (zmScheme_)
{ {
case ZmScheme::qedTL: case ZmScheme::qedTL:
{ {
std::vector<int> zm(nd,0); std::vector<int> zm(grid_->Nd(), 0);
TComplex Tzero = Complex(0.0,0.0); ScalarSite z = ScalarComplex(0., 0.);
pokeSite(Tzero, out, zm);
pokeSite(z, out, zm);
break; break;
} }
case ZmScheme::qedL: case ZmScheme::qedL:
{ {
LatticeInteger spNrm(grid), coor(grid); LatticeInteger spNrm(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);
// IR improvement makeSpatialNorm(spNrm);
out = where(spNrm == Integer(0), 0.*out, out);
for(int i = 0; i < improvement_.size(); i++) for(int i = 0; i < improvement_.size(); i++)
{ {
Real f = sqrt(improvement_[i]+1); Real f = sqrt(improvement_[i] + 1);
out = where(spNrm == Integer(i+1), f*out, out); out = where(spNrm == Integer(i + 1), f*out, out);
} }
} }
default: default:
assert(0);
break; break;
} }
} }
template<class Gimpl> template<class GImpl>
void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in, void Photon<GImpl>::transverseProjectSpatial(GaugeField &out)
GaugeField &out)
{ {
GridBase *grid = out._grid; const unsigned int nd = grid_->Nd();
LatticeComplex momProp(grid); GaugeLinkField invKHat(grid_), cst(grid_);
LatticeInteger spNrm(grid_);
switch (zmScheme_) std::vector<GaugeLinkField> khat, a(nd, grid_), aProj(nd, grid_);
invKHat = zero;
makeSpatialNorm(spNrm);
makeKHat(khat);
for (unsigned int mu = 0; mu < nd; ++mu)
{ {
case ZmScheme::qedTL: a[mu] = peekLorentz(out, mu);
case ZmScheme::qedL: if (mu < nd - 1)
{ {
invKHatSquared(momProp); invKHat += khat[mu]*conjugate(khat[mu]);
zmSub(momProp);
break;
} }
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); aProj[mu] -= invKHat*khat[mu]*conjugate(khat[nu])*a[nu];
break;
} }
pokeLorentz(out, aProj[mu], mu);
}
}
template<class GImpl>
void Photon<GImpl>::gaugeTransform(GaugeField &out)
{
switch (gauge_)
{
case Gauge::feynman:
break;
case Gauge::coulomb:
transverseProjectSpatial(out);
break;
case Gauge::landau:
assert(0);
break;
default: default:
assert(0);
break; break;
} }
}
template<class GImpl>
void Photon<GImpl>::MomentumSpacePropagator(const GaugeField &in,
GaugeField &out)
{
LatticeComplex momProp(grid_);
makeInvKHatSquared(momProp);
zmSub(momProp);
out = in*momProp; out = in*momProp;
} }
template<class Gimpl> template<class GImpl>
void Photon<Gimpl>::StochasticWeight(GaugeLinkField &weight) void Photon<GImpl>::StochasticWeight(GaugeLinkField &weight)
{ {
auto *grid = dynamic_cast<GridCartesian *>(weight._grid); const unsigned int nd = grid_->Nd();
const unsigned int nd = grid->_ndimension; std::vector<int> l = grid_->FullDimensions();
std::vector<int> latt_size = grid->_fdimensions; Integer vol = 1;
switch (zmScheme_) for(unsigned int mu = 0; mu < nd; mu++)
{ {
case ZmScheme::qedTL: vol = vol*l[mu];
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;
} }
makeInvKHatSquared(weight);
weight = sqrt(vol)*sqrt(weight);
zmSub(weight);
} }
template<class Gimpl> template<class GImpl>
void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng) void Photon<GImpl>::StochasticField(GaugeField &out, GridParallelRNG &rng)
{ {
auto *grid = dynamic_cast<GridCartesian *>(out._grid); GaugeLinkField weight(grid_);
GaugeLinkField weight(grid);
StochasticWeight(weight); StochasticWeight(weight);
StochasticField(out, rng, weight); StochasticField(out, rng, weight);
} }
template<class Gimpl> template<class GImpl>
void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng, void Photon<GImpl>::StochasticField(GaugeField &out, GridParallelRNG &rng,
const GaugeLinkField &weight) const GaugeLinkField &weight)
{ {
auto *grid = dynamic_cast<GridCartesian *>(out._grid); const unsigned int nd = grid_->Nd();
const unsigned int nd = grid->_ndimension; GaugeLinkField r(grid_);
GaugeLinkField r(grid); GaugeField aTilde(grid_);
GaugeField aTilde(grid); FFT fft(dynamic_cast<GridCartesian *>(grid_));
FFT fft(grid);
switch (zmScheme_) for(unsigned int mu = 0; mu < nd; mu++)
{ {
case ZmScheme::qedTL: gaussian(rng, r);
case ZmScheme::qedL: r = weight*r;
{ pokeLorentz(aTilde, r, mu);
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;
} }
gaugeTransform(aTilde);
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>
void Photon<Gimpl>::UnitField(GaugeField &out) void Photon<GImpl>::UnitField(GaugeField &out)
{ {
auto *grid = dynamic_cast<GridCartesian *>(out._grid); const unsigned int nd = grid_->Nd();
const unsigned int nd = grid->_ndimension; GaugeLinkField r(grid_);
GaugeLinkField r(grid);
r = Complex(1.0,0.0); r = ScalarComplex(1., 0.);
for(unsigned int mu = 0; mu < nd; mu++)
for(int mu = 0; mu < nd; mu++)
{ {
pokeLorentz(out, r, mu); pokeLorentz(out, r, mu);
} }
out = real(out); out = real(out);
} }
// template<class Gimpl>
// void Photon<Gimpl>::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<class Gimpl>
// 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

2
Hadrons/Modules/MGauge/StochEm.cc
View File

@ -70,7 +70,7 @@ void TStochEm::execute(void)
LOG(Message) << "Generating stochastic EM potential..." << std::endl; LOG(Message) << "Generating stochastic EM potential..." << std::endl;
std::vector<Real> improvements = strToVec<Real>(par().improvement); std::vector<Real> improvements = strToVec<Real>(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 &a = envGet(EmField, getName());
auto &w = envGet(EmComp, "_" + getName() + "_weight"); auto &w = envGet(EmComp, "_" + getName() + "_weight");

3
Hadrons/Modules/MGauge/StochEm.hpp
View File

@ -47,8 +47,7 @@ public:
GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar, GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar,
PhotonR::Gauge, gauge, PhotonR::Gauge, gauge,
PhotonR::ZmScheme, zmScheme, PhotonR::ZmScheme, zmScheme,
std::string, improvement, std::string, improvement);
Real, G0_qedInf);
}; };
class TStochEm: public Module<StochEmPar> class TStochEm: public Module<StochEmPar>

2
Hadrons/Modules/MGauge/UnitEm.cc
View File

@ -62,7 +62,7 @@ void TUnitEm::setup(void)
// execution /////////////////////////////////////////////////////////////////// // execution ///////////////////////////////////////////////////////////////////
void TUnitEm::execute(void) 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()); auto &a = envGet(EmField, getName());
LOG(Message) << "Generating unit EM potential..." << std::endl; LOG(Message) << "Generating unit EM potential..." << std::endl;
photon.UnitField(a); photon.UnitField(a);