diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index ed6254f4..734ff0fc 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include #include #include #include diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index ee8d8d56..0fb23cc4 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp +Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp Copyright (C) 2015-2018 Author: Antonin Portelli -Author: Lanny91 +Author: Andrew Lawson +Author: Vera Guelpers 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 @@ -34,13 +35,17 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +#include + BEGIN_HADRONS_NAMESPACE /* - Sequential source + Sequential source with insertion of conserved current. + Additionally optional insertion of a photon field A_\mu(x). ----------------------------- - * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) + * src_x = sum_{mu=mu_min}^{mu_max} + q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) (* A_\mu(x)) * options: - q: input propagator (string) @@ -48,8 +53,10 @@ BEGIN_HADRONS_NAMESPACE - tA: begin timeslice (integer) - tB: end timesilce (integer) - curr_type: type of conserved current to insert (Current) - - mu: Lorentz index of current to insert (integer) + - mu_min: begin Lorentz Index (integer) + - mu_max: end Lorentz Index (integer) - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + - photon: optional photon field (string) */ @@ -67,8 +74,10 @@ public: unsigned int, tA, unsigned int, tB, Current, curr_type, - unsigned int, mu, - std::string, mom); + unsigned int, mu_min, + unsigned int, mu_max, + std::string, mom, + std::string, photon); }; template @@ -76,6 +85,8 @@ class TSeqConserved: public Module { public: FERM_TYPE_ALIASES(FImpl,); +public: + typedef PhotonR::GaugeField EmField; public: // constructor TSeqConserved(const std::string name); @@ -93,6 +104,7 @@ protected: MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); + /****************************************************************************** * TSeqConserved implementation * ******************************************************************************/ @@ -107,7 +119,8 @@ template std::vector TSeqConserved::getInput(void) { std::vector in = {par().q, par().action}; - + if (!par().photon.empty()) in.push_back(par().photon); + return in; } @@ -116,7 +129,7 @@ std::vector TSeqConserved::getOutput(void) { std::vector out = {getName()}; - return out; + return out; } // setup /////////////////////////////////////////////////////////////////////// @@ -125,6 +138,10 @@ void TSeqConserved::setup(void) { auto Ls_ = env().getObjectLs(par().action); envCreateLat(PropagatorField, getName(), Ls_); + envTmpLat(PropagatorField, "src_tmp"); + envTmpLat(LatticeComplex, "mom_phase"); + envTmpLat(LatticeComplex, "coor"); + envTmpLat(LatticeComplex, "latt_compl"); } // execution /////////////////////////////////////////////////////////////////// @@ -134,27 +151,71 @@ void TSeqConserved::execute(void) if (par().tA == par().tB) { LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion (mu = " - << par().mu << ") at " << "t = " << par().tA << std::endl; + << par().curr_type << " current at " + << "t = " << par().tA << std::endl; } else { LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion (mu = " - << par().mu << ") for " << par().tA << " <= t <= " + << par().curr_type << " current for " + << par().tA << " <= t <= " << par().tB << std::endl; } auto &src = envGet(PropagatorField, getName()); + envGetTmp(PropagatorField, src_tmp); + src_tmp = src; auto &q = envGet(PropagatorField, par().q); auto &mat = envGet(FMat, par().action); + envGetTmp(LatticeComplex, mom_phase); + envGetTmp(LatticeComplex, coor); + envGetTmp(LatticeComplex, latt_compl); + src = zero; + + //exp(ipx) std::vector mom = strToVec(par().mom); - mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, - mom, par().tA, par().tB); + mom_phase = zero; + Complex i(0.0,1.0); + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor; + } + mom_phase = exp((Real)(2*M_PI)*i*mom_phase); + LOG(Message) << "Inserting momentum " << mom << std::endl; + + + + if (!par().photon.empty()) + { + LOG(Message) << "Inserting the stochastic photon field " << par().photon << std::endl; + } + + for(unsigned int mu=par().mu_min;mu<=par().mu_max;mu++) + { + if (!par().photon.empty()) + { + //Get the stochastic photon field, if required + auto &stoch_photon = envGet(EmField, par().photon); + latt_compl = PeekIndex(stoch_photon, mu) * mom_phase; + } + else + { + latt_compl = mom_phase; + } + + mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, + par().tA, par().tB, latt_compl); + src += src_tmp; + + } + + } + END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_SeqConserved_hpp_ +#endif // Hadrons_MSource_SeqConserved_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp deleted file mode 100644 index 243d05f2..00000000 --- a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp +++ /dev/null @@ -1,169 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp - -Copyright (C) 2015-2018 - -Author: Antonin Portelli -Author: Andrew Lawson -Author: Vera Guelpers - -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 Hadrons_MSource_SeqConservedSummed_hpp_ -#define Hadrons_MSource_SeqConservedSummed_hpp_ - -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/* - - Sequential source summed over the Lorentz index of current - ----------------------------- - * src_x = sum_mu q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) - - * options: - - q: input propagator (string) - - action: fermion action used for propagator q (string) - - tA: begin timeslice (integer) - - tB: end timesilce (integer) - - curr_type: type of conserved current to insert (Current) - - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") - - */ - -/****************************************************************************** - * SeqConservedSummed * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MSource) - -class SeqConservedSummedPar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedSummedPar, - std::string, q, - std::string, action, - unsigned int, tA, - unsigned int, tB, - Current, curr_type, - std::string, mom); -}; - -template -class TSeqConservedSummed: public Module -{ -public: - FERM_TYPE_ALIASES(FImpl,); -public: - // constructor - TSeqConservedSummed(const std::string name); - // destructor - virtual ~TSeqConservedSummed(void) = default; - // dependency relation - virtual std::vector getInput(void); - virtual std::vector getOutput(void); -protected: - // setup - virtual void setup(void); - // execution - virtual void execute(void); -}; - -MODULE_REGISTER_NS(SeqConservedSummed, TSeqConservedSummed, MSource); - -/****************************************************************************** - * TSeqConservedSummed implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TSeqConservedSummed::TSeqConservedSummed(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TSeqConservedSummed::getInput(void) -{ - std::vector in = {par().q, par().action}; - - return in; -} - -template -std::vector TSeqConservedSummed::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TSeqConservedSummed::setup(void) -{ - auto Ls_ = env().getObjectLs(par().action); - envCreateLat(PropagatorField, getName(), Ls_); - envTmpLat(PropagatorField, "src_tmp"); -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TSeqConservedSummed::execute(void) -{ - if (par().tA == par().tB) - { - LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion summed over mu at " - << "t = " << par().tA << std::endl; - } - else - { - LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion summed over mu for " - << par().tA << " <= t <= " - << par().tB << std::endl; - } - auto &src = envGet(PropagatorField, getName()); - envGetTmp(PropagatorField, src_tmp); - src_tmp = src; - auto &q = envGet(PropagatorField, par().q); - auto &mat = envGet(FMat, par().action); - - std::vector mom = strToVec(par().mom); - src = zero; - for(int mu=0;mu<=3;mu++) - { - mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, - mom, par().tA, par().tB); - src += src_tmp; - - } - -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_SeqConservedSummed_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index d005caec..63b674d4 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -20,7 +20,6 @@ modules_hpp =\ Modules/MContraction/WardIdentity.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MFermion/GaugeProp.hpp \ - Modules/MSource/SeqConservedSummed.hpp \ Modules/MSource/SeqGamma.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/Wall.hpp \ diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 1d395d53..94e065cc 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -125,9 +125,9 @@ namespace Grid { PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax)=0; + unsigned int tmax, + Lattice> &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 5005b1a2..3a296e52 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -407,17 +407,19 @@ void ImprovedStaggeredFermion::ContractConservedCurrent(PropagatorField &q } 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) +void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx) { assert(0); + } + FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 9d5270c6..015eb6bb 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -166,13 +166,13 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 7308ff30..3640a26b 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -419,15 +419,16 @@ void ImprovedStaggeredFermion5D::ContractConservedCurrent(PropagatorField } 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) +void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx) { assert(0); + } FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 22fb9e7d..1c211938 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -178,13 +178,13 @@ namespace QCD { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 1a020e8a..560c371d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -381,40 +381,100 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, } } +//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); +// PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); +// unsigned int tshift = (mu == Tp) ? 1 : 0; +// unsigned int LLt = GridDefaultLatt()[Tp]; +// +// // 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); +// +// q_out = zero; +// LatticeInteger coords(_grid); +// LatticeCoordinate(coords, Tp); +// +// // Need q(x + mu) and q(x - mu). +// tmp = Cshift(q_in, mu, 1); +// tmpFwd = tmp*ph; +// tmp = ph*q_in; +// tmpBwd = Cshift(tmp, mu, -1); +// +// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) +// { +// // Compute the sequential conserved current insertion only if our simd +// // object contains a timeslice we need. +// vInteger t_mask = ((coords._odata[sU] >= tmin) && +// (coords._odata[sU] <= tmax)); +// Integer timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], +// q_out._odata[sU], +// Umu, sU, mu, t_mask); +// } +// +// // Repeat for backward direction. +// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && +// (coords._odata[sU] <= (tmax + tshift))); +// +// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) +// unsigned int t0 = 0; +// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); +// +// timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], +// q_out._odata[sU], +// Umu, sU, mu, t_mask); +// } +// } +//} + 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) + unsigned int tmax, + Lattice> &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); - Lattice> ph(_grid), coor(_grid); Complex i(0.0,1.0); PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLt = GridDefaultLatt()[Tp]; - // 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); - q_out = zero; LatticeInteger coords(_grid); LatticeCoordinate(coords, Tp); // Need q(x + mu) and q(x - mu). tmp = Cshift(q_in, mu, 1); - tmpFwd = tmp*ph; - tmp = ph*q_in; + tmpFwd = tmp*lattice_cmplx; + tmp = lattice_cmplx*q_in; tmpBwd = Cshift(tmp, mu, -1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) @@ -449,6 +509,8 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, Umu, sU, mu, t_mask); } } + + } FermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index feba40ed..db32529d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -155,13 +155,13 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { 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); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index bc703187..92280c13 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,18 +779,110 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, } +//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(FermionGrid()), coor(FermionGrid()); +// PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), +// tmp(FermionGrid()); +// Complex i(0.0, 1.0); +// unsigned int tshift = (mu == Tp) ? 1 : 0; +// unsigned int LLs = q_in._grid->_rdimensions[0]; +// unsigned int LLt = GridDefaultLatt()[Tp]; +// +// // Momentum projection. +// ph = zero; +// for(unsigned int nu = 0; nu < Nd - 1; nu++) +// { +// // Shift coordinate lattice index by 1 to account for 5th dimension. +// LatticeCoordinate(coor, nu + 1); +// ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); +// } +// ph = exp((Real)(2*M_PI)*i*ph); +// +// q_out = zero; +// LatticeInteger coords(_FourDimGrid); +// LatticeCoordinate(coords, Tp); +// +// +// // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu +// // by one. +// tmp = Cshift(q_in, mu + 1, 1); +// tmpFwd = tmp*ph; +// tmp = ph*q_in; +// tmpBwd = Cshift(tmp, mu + 1, -1); +// +// parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) +// { +// // Compute the sequential conserved current insertion only if our simd +// // object contains a timeslice we need. +// vInteger t_mask = ((coords._odata[sU] >= tmin) && +// (coords._odata[sU] <= tmax)); +// Integer timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// unsigned int sF = sU * LLs; +// for (unsigned int s = 0; s < LLs; ++s) +// { +// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); +// bool tadpole_sign = (curr_type == Current::Tadpole); +// bool switch_sgn = tadpole_sign || axial_sign; +// +// Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], +// q_out._odata[sF], Umu, sU, +// mu, t_mask, switch_sgn); +// ++sF; +// } +// } +// +// // Repeat for backward direction. +// t_mask = ((coords._odata[sU] >= (tmin + tshift)) && +// (coords._odata[sU] <= (tmax + tshift))); +// +// //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3) +// unsigned int t0 = 0; +// if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 )); +// +// timeSlices = Reduce(t_mask); +// +// if (timeSlices > 0) +// { +// unsigned int sF = sU * LLs; +// for (unsigned int s = 0; s < LLs; ++s) +// { +// bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); +// Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], +// q_out._odata[sF], Umu, sU, +// mu, t_mask, axial_sign); +// ++sF; +// } +// } +// } +//} + + + + 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) + unsigned int tmax, + Lattice> &lattice_cmplx) { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); - Lattice> ph(FermionGrid()), coor(FermionGrid()); PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), tmp(FermionGrid()); Complex i(0.0, 1.0); @@ -798,25 +890,26 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, unsigned int LLs = q_in._grid->_rdimensions[0]; unsigned int LLt = GridDefaultLatt()[Tp]; - // Momentum projection. - ph = zero; - for(unsigned int nu = 0; nu < Nd - 1; nu++) - { - // Shift coordinate lattice index by 1 to account for 5th dimension. - LatticeCoordinate(coor, nu + 1); - ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); - } - ph = exp((Real)(2*M_PI)*i*ph); - q_out = zero; LatticeInteger coords(_FourDimGrid); LatticeCoordinate(coords, Tp); + + //QED: photon field is 4dim, but need a 5dim object to multiply to + // DWF PropagatorField + Lattice> lattice_cmplx_5d(FermionGrid()); + for (unsigned int s = 0; s < LLs; ++s) + { + InsertSlice(lattice_cmplx,lattice_cmplx_5d, s, 0); + } + + + // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu // by one. tmp = Cshift(q_in, mu + 1, 1); - tmpFwd = tmp*ph; - tmp = ph*q_in; + tmpFwd = tmp*lattice_cmplx_5d; + tmp = lattice_cmplx_5d*q_in; tmpBwd = Cshift(tmp, mu + 1, -1); parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) @@ -868,6 +961,10 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, } } + + + + FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d66f4a1d..6a7c6c7d 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -222,13 +222,13 @@ namespace QCD { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, + void SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, - Current curr_type, + Current curr_type, unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + unsigned int tmin, + unsigned int tmax, + Lattice> &lattice_cmplx); }; }} diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 0265f5a6..6188a580 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/tests/hadrons/Test_hadrons.hpp @@ -263,7 +263,8 @@ inline void makeConservedSequentialSource(Application &application, seqPar.tA = tS; seqPar.tB = tS; seqPar.curr_type = curr; - seqPar.mu = mu; + seqPar.mu_min = mu; + seqPar.mu_min = mu; seqPar.mom = mom; application.createModule(srcName, seqPar); }