diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 6e6144da..c66c7b13 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -489,6 +489,14 @@ namespace QCD { return traceIndex(lhs); } + ////////////////////////////////////////// + // Current types + ////////////////////////////////////////// + GRID_SERIALIZABLE_ENUM(Current, undef, + Vector, 0, + Axial, 1, + Tadpole, 2); + } //namespace QCD } // Grid diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 676a0e83..144b70f6 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -112,6 +112,21 @@ namespace Grid { /////////////////////////////////////////////// virtual void ImportGauge(const GaugeField & _U)=0; + ////////////////////////////////////////////////////////////////////// + // Conserved currents, either contract at sink or insert sequentially. + ////////////////////////////////////////////////////////////////////// + virtual void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu)=0; + virtual void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 2ba4f4af..ef8c79bd 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -395,6 +395,31 @@ void ImprovedStaggeredFermion::DhopInternal(StencilImpl &st, LebesgueOrder } }; +//////////////////////////////////////////////////////// +// Conserved current - not yet implemented. +//////////////////////////////////////////////////////// +template +void ImprovedStaggeredFermion::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + assert(0); +} + FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 7d1f2996..9d5270c6 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -157,6 +157,22 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 61a3c559..293077f7 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -345,6 +345,30 @@ void ImprovedStaggeredFermion5D::MooeeInvDag(const FermionField &in, MooeeInv(in, out); } +//////////////////////////////////////////////////////// +// Conserved current - not yet implemented. +//////////////////////////////////////////////////////// +template +void ImprovedStaggeredFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + assert(0); +} + +template +void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + assert(0); +} FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); FermOpStaggeredVec5dTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 4961da49..1c540892 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -160,6 +160,21 @@ namespace QCD { // Comms buffer std::vector > comm_buf; + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 32083d5e..839f5215 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -347,6 +347,53 @@ void WilsonFermion::DhopInternal(StencilImpl &st, LebesgueOrder &lo, } }; +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * to make a conserved current sink or inserting the conserved current + * sequentially. + ******************************************************************************/ +template +void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + Gamma g5(Gamma::Algebra::Gamma5); + conformable(_grid, q_in_1._grid); + conformable(_grid, q_in_2._grid); + conformable(_grid, q_out._grid); + Kernels::ContractConservedCurrentInternal(q_in_1, q_in_2, q_out, + Umu, curr_type, mu); +} + +template +void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + conformable(_grid, q_in._grid); + conformable(_grid, q_out._grid); + Lattice> ph(_grid), coor(_grid); + Complex i(0.0,1.0); + + // Momentum projection + ph = zero; + for(unsigned int mu = 0; mu < Nd - 1; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + + Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, mu, ph, + tmin, tmax); +} + FermOpTemplateInstantiate(WilsonFermion); AdjointFermOpTemplateInstantiate(WilsonFermion); TwoIndexFermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 933be732..feba40ed 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -146,6 +146,22 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 88bc425a..d0d3d055 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -679,6 +679,83 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe } +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * to make a conserved current sink or inserting the conserved current + * sequentially. + ******************************************************************************/ +template +void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu) +{ + conformable(q_in_1._grid, FermionGrid()); + conformable(q_in_1._grid, q_in_2._grid); + conformable(_FourDimGrid, q_out._grid); + + PropagatorField q1_s(_FourDimGrid); + PropagatorField q2_s(_FourDimGrid); + PropagatorField tmp(_FourDimGrid); + + // Contract across 5th dimension. + q_out = zero; + for (int s = 0; s < Ls; ++s) + { + ExtractSlice(q1_s, q_in_1, 0, s); + ExtractSlice(q2_s, q_in_2, 0, Ls - s - 1); + Kernels::ContractConservedCurrentInternal(q1_s, q2_s, tmp, Umu, curr_type, mu); + + // Axial current sign + Real G_s = (curr_type == Current::Axial) ? ((s < Ls/2) ? -1. : 1.) : 1.; + q_out += G_s*tmp; + } +} + + +template +void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax) +{ + conformable(q_in._grid, FermionGrid()); + conformable(q_in._grid, q_out._grid); + Lattice> ph(_FourDimGrid), coor(_FourDimGrid); + Complex i(0.0, 1.0); + + // Momentum projection + ph = zero; + for(unsigned int nu = 0; nu < Nd - 1; nu++) + { + LatticeCoordinate(coor, nu); + ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + + // Sequential insertion + Kernels::SeqConservedCurrentInternal(q_in, q_out, Umu, curr_type, + mu, ph, tmin, tmax); + + // Axial current sign. + if (curr_type == Current::Axial) + { + SitePropagator result; + parallel_for(int sU = 0; sU < Umu._grid->oSites(); sU++) + { + int sF = sU * Ls; + for (int s = 0; s < Ls/2; s++) + { + vstream(q_out._odata[sF], -q_out._odata[sF]); + sF++; + } + } + } +} FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index e87e927e..d66f4a1d 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -214,6 +214,21 @@ namespace QCD { // Comms buffer std::vector > comm_buf; + /////////////////////////////////////////////////////////////// + // Conserved current utilities + /////////////////////////////////////////////////////////////// + void ContractConservedCurrent(PropagatorField &q_in_1, + PropagatorField &q_in_2, + PropagatorField &q_out, + Current curr_type, + unsigned int mu); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + std::vector mom, + unsigned int tmin, + unsigned int tmax); }; }} diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 6e72e089..fbf8dc00 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -567,6 +567,95 @@ void WilsonKernels::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal vstream(out._odata[sF], result); } +/******************************************************************************* + * Conserved current utilities for Wilson fermions, for contracting propagators + * to make a conserved current sink or inserting the conserved current + * sequentially. Common to both 4D and 5D. + ******************************************************************************/ +#define WilsonCurrentFwd(expr, mu) (0.5*(Gamma::gmu[mu]*expr - expr)) +#define WilsonCurrentBwd(expr, mu) (0.5*(Gamma::gmu[mu]*expr + expr)) + +template +void WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu) +{ + Gamma g5(Gamma::Algebra::Gamma5); + PropagatorField tmp(q_out._grid); + GaugeLinkField Umu(U._grid); + Umu = PeekIndex(U, mu); + + tmp = this->CovShiftForward(Umu, mu, q_in_1); + q_out = (g5*adj(q_in_2)*g5)*WilsonCurrentFwd(tmp, mu); + + tmp = adj(Umu)*q_in_1; + q_out += (g5*adj(this->CovShiftForward(Umu, mu, q_in_2))*g5)*WilsonCurrentBwd(q_in_1, mu); +} + + +template +void WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu, + Lattice> &ph, + unsigned int tmin, + unsigned int tmax) +{ + int tshift = (mu == Nd - 1) ? 1 : 0; + Real G_T = (curr_type == Current::Tadpole) ? -1. : 1.; + PropagatorField tmp(q_in._grid); + GaugeLinkField Umu(U._grid); + Umu = PeekIndex(U, mu); + Lattice> t(q_in._grid); + + tmp = this->CovShiftForward(Umu, mu, q_in)*ph; + where((t >= tmin) and (t <= tmax), tmp, 0.*tmp); + q_out = G_T*WilsonCurrentFwd(tmp, mu); + + tmp = q_in*ph; + tmp = this->CovShiftBackward(Umu, mu, tmp); + where((t >= tmin + tshift) and (t <= tmax + tshift), tmp, 0.*tmp); + q_out += WilsonCurrentBwd(tmp, mu); +} + + +// GParity, (Z)DomainWallVec5D -> require special implementation +#define NO_CURR(Impl) \ +template <> void \ +WilsonKernels::ContractConservedCurrentInternal(const PropagatorField &q_in_1, \ + const PropagatorField &q_in_2, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + Current curr_type, \ + unsigned int mu) \ +{ \ + assert(0); \ +} \ +template <> void \ +WilsonKernels::SeqConservedCurrentInternal(const PropagatorField &q_in, \ + PropagatorField &q_out, \ + DoubledGaugeField &U, \ + Current curr_type, \ + unsigned int mu, \ + Lattice> &ph, \ + unsigned int tmin, \ + unsigned int tmax) \ +{ \ + assert(0); \ +} + +NO_CURR(GparityWilsonImplF); +NO_CURR(GparityWilsonImplD); +NO_CURR(DomainWallVec5dImplF); +NO_CURR(DomainWallVec5dImplD); +NO_CURR(ZDomainWallVec5dImplF); +NO_CURR(ZDomainWallVec5dImplD); + FermOpTemplateInstantiate(WilsonKernels); AdjointFermOpTemplateInstantiate(WilsonKernels); TwoIndexFermOpTemplateInstantiate(WilsonKernels); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 20ee87f2..34820274 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -166,6 +166,24 @@ public: void DhopDir(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf, int sF, int sU, const FermionField &in, FermionField &out, int dirdisp, int gamma); + ////////////////////////////////////////////////////////////////////////////// + // Utilities for inserting Wilson conserved current. + ////////////////////////////////////////////////////////////////////////////// + void ContractConservedCurrentInternal(const PropagatorField &q_in_1, + const PropagatorField &q_in_2, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu); + void SeqConservedCurrentInternal(const PropagatorField &q_in, + PropagatorField &q_out, + DoubledGaugeField &U, + Current curr_type, + unsigned int mu, + Lattice> &ph, + unsigned int tmin, + unsigned int tmax); + private: // Specialised variants void GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,