From 55e39df30f1953e0d54b79d2b1c6061093468e76 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 22 Dec 2017 11:36:31 +0000 Subject: [PATCH 1/8] tadpole insertion for DWF --- lib/qcd/action/fermion/WilsonFermion5D.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 393ee7f3..bc703187 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -833,9 +833,12 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, 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, axial_sign); + mu, t_mask, switch_sgn); ++sF; } } From 935cd1e173e6f50f5185efb30e6b9c3b92790c3d Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 22 Dec 2017 11:38:45 +0000 Subject: [PATCH 2/8] conserved current insertion summed over Lorentzindex --- extras/Hadrons/Modules.hpp | 67 ++----- .../Modules/MSource/SeqConservedSummed.hpp | 165 ++++++++++++++++++ extras/Hadrons/modules.inc | 46 ++--- 3 files changed, 208 insertions(+), 70 deletions(-) create mode 100644 extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index e1f06f32..a52b9a6e 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,60 +1,31 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules.hpp - -Copyright (C) 2015 -Copyright (C) 2016 -Copyright (C) 2017 - -Author: Antonin Portelli - -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 */ - -#include -#include -#include -#include -#include -#include +#include #include -#include +#include +#include +#include #include -#include +#include +#include #include -#include -#include -#include -#include -#include +#include +#include +#include #include -#include #include +#include #include #include #include -#include +#include +#include +#include #include #include #include -#include #include -#include -#include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp new file mode 100644 index 00000000..2296701a --- /dev/null +++ b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp @@ -0,0 +1,165 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp + +Copyright (C) 2017 + +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); + // 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); + env().template registerLattice(getName(), Ls_); +} + +// 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; + } + PropagatorField &src = *env().template createLattice(getName()); + auto src_buf = src; + PropagatorField &q = *env().template getObject(par().q); + FMat &mat = *(env().template getObject(par().action)); + + std::vector mom = strToVec(par().mom); + src = zero; + for(int mu=0;mu<=3;mu++) + { + mat.SeqConservedCurrent(q, src_buf, par().curr_type, mu, + mom, par().tA, par().tB); + src += src_buf; + + } + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_SeqConservedSummed_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index fbbb2eb9..1ad4750a 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,42 +1,44 @@ modules_cc =\ Modules/MContraction/WeakHamiltonianEye.cc \ - Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ - Modules/MGauge/Load.cc \ + Modules/MContraction/WeakHamiltonianNonEye.cc \ + Modules/MScalar/FreeProp.cc \ + Modules/MScalar/ChargedProp.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Unit.cc \ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/FreeProp.cc + Modules/MGauge/Load.cc modules_hpp =\ - Modules/MAction/DWF.hpp \ - Modules/MAction/Wilson.hpp \ - Modules/MContraction/Baryon.hpp \ - Modules/MContraction/DiscLoop.hpp \ - Modules/MContraction/Gamma3pt.hpp \ - Modules/MContraction/Meson.hpp \ + Modules/MSolver/RBPrecCG.hpp \ Modules/MContraction/WardIdentity.hpp \ - Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MContraction/Gamma3pt.hpp \ + Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ - Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MContraction/Baryon.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ Modules/MContraction/WeakNeutral4ptDisc.hpp \ - Modules/MFermion/GaugeProp.hpp \ - Modules/MGauge/Load.hpp \ - Modules/MGauge/Random.hpp \ - Modules/MGauge/StochEm.hpp \ - Modules/MGauge/Unit.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MUtilities/TestSeqGamma.hpp \ + Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ - Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ + Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/Scalar.hpp \ Modules/MSink/Point.hpp \ Modules/MSink/Smear.hpp \ - Modules/MSolver/RBPrecCG.hpp \ + Modules/MFermion/GaugeProp.hpp \ + Modules/MSource/SeqConservedSummed.hpp \ + Modules/MSource/Wall.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqConserved.hpp \ Modules/MSource/SeqGamma.hpp \ - Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ - Modules/MUtilities/TestSeqConserved.hpp \ - Modules/MUtilities/TestSeqGamma.hpp + Modules/MGauge/Load.hpp \ + Modules/MGauge/StochEm.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MAction/DWF.hpp \ + Modules/MAction/Wilson.hpp + From 389731d373e5f348c80a29520a4328414470c0a2 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 23 Jan 2018 10:11:33 +0000 Subject: [PATCH 3/8] changed SeqConservedSummed.hpp to work with new hadrons interface --- .../Modules/MSource/SeqConservedSummed.hpp | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp index 2296701a..243d05f2 100644 --- a/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp @@ -4,8 +4,9 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp -Copyright (C) 2017 +Copyright (C) 2015-2018 +Author: Antonin Portelli Author: Andrew Lawson Author: Vera Guelpers @@ -82,6 +83,7 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); +protected: // setup virtual void setup(void); // execution @@ -121,7 +123,8 @@ template void TSeqConservedSummed::setup(void) { auto Ls_ = env().getObjectLs(par().action); - env().template registerLattice(getName(), Ls_); + envCreateLat(PropagatorField, getName(), Ls_); + envTmpLat(PropagatorField, "src_tmp"); } // execution /////////////////////////////////////////////////////////////////// @@ -141,18 +144,19 @@ void TSeqConservedSummed::execute(void) << par().tA << " <= t <= " << par().tB << std::endl; } - PropagatorField &src = *env().template createLattice(getName()); - auto src_buf = src; - PropagatorField &q = *env().template getObject(par().q); - FMat &mat = *(env().template getObject(par().action)); + 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_buf, par().curr_type, mu, + mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, mom, par().tA, par().tB); - src += src_buf; + src += src_tmp; } From b6fe03eb26a8725d4f52d3215104ea9743c071b0 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 2 Feb 2018 15:29:38 +0000 Subject: [PATCH 4/8] BugFix: Now the stochatic EM potential weight is generated when calling for the first time --- extras/Hadrons/Modules/MGauge/StochEm.cc | 6 ++++-- extras/Hadrons/Modules/MGauge/StochEm.hpp | 2 ++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc index 21b7f626..cfbd28ef 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.cc +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -57,9 +57,11 @@ std::vector TStochEm::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TStochEm::setup(void) { + create_weight = false; if (!env().hasCreatedObject("_" + getName() + "_weight")) { envCacheLat(EmComp, "_" + getName() + "_weight"); + create_weight = true; } envCreateLat(EmField, getName()); } @@ -67,13 +69,13 @@ void TStochEm::setup(void) // execution /////////////////////////////////////////////////////////////////// void TStochEm::execute(void) { - LOG(Message) << "Generating stochatic EM potential..." << std::endl; + LOG(Message) << "Generating stochastic EM potential..." << std::endl; PhotonR photon(par().gauge, par().zmScheme); auto &a = envGet(EmField, getName()); auto &w = envGet(EmComp, "_" + getName() + "_weight"); - if (!env().hasCreatedObject("_" + getName() + "_weight")) + if (create_weight) { LOG(Message) << "Caching stochatic EM potential weight (gauge: " << par().gauge << ", zero-mode scheme: " diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 87b70880..05162443 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -60,6 +60,8 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); +private: + bool create_weight; protected: // setup virtual void setup(void); From 49a0ae73eb097fc9e7de02897bcde6ff6a3db170 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Mon, 12 Feb 2018 09:36:08 +0000 Subject: [PATCH 5/8] Insertion of photon field in seqential conserved current --- extras/Hadrons/Modules.hpp | 1 - .../Hadrons/Modules/MSource/SeqConserved.hpp | 93 ++++++++-- .../Modules/MSource/SeqConservedSummed.hpp | 169 ------------------ extras/Hadrons/modules.inc | 1 - lib/qcd/action/fermion/FermionOperator.h | 4 +- .../fermion/ImprovedStaggeredFermion.cc | 16 +- .../action/fermion/ImprovedStaggeredFermion.h | 12 +- .../fermion/ImprovedStaggeredFermion5D.cc | 15 +- .../fermion/ImprovedStaggeredFermion5D.h | 12 +- lib/qcd/action/fermion/WilsonFermion.cc | 90 ++++++++-- lib/qcd/action/fermion/WilsonFermion.h | 14 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 127 +++++++++++-- lib/qcd/action/fermion/WilsonFermion5D.h | 10 +- tests/hadrons/Test_hadrons.hpp | 3 +- 14 files changed, 310 insertions(+), 257 deletions(-) delete mode 100644 extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp 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); } From c6e1f645737d6785d513246db47323afaa5850e2 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Tue, 13 Feb 2018 09:30:23 +0000 Subject: [PATCH 6/8] Test for QED --- tests/hadrons/Test_QED.cc | 260 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 tests/hadrons/Test_QED.cc diff --git a/tests/hadrons/Test_QED.cc b/tests/hadrons/Test_QED.cc new file mode 100644 index 00000000..44dd0578 --- /dev/null +++ b/tests/hadrons/Test_QED.cc @@ -0,0 +1,260 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015-2018 + + Author: Antonin Portelli + 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. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.2}; //{.01, .04, .2 , .25 , .3 }; + + unsigned int nt = GridDefaultLatt()[Tp]; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + // pt source + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + //stochastic photon field + MGauge::StochEm::Par photonPar; + photonPar.gauge = PhotonR::Gauge::feynman; + photonPar.zmScheme = PhotonR::ZmScheme::qedL; + application.createModule("ph_field", photonPar); + + + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 8; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + application.createModule("DWF_" + flavour[i], actionPar); + + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], + quarkPar); + + + //seq sources with tadpole insertion + MSource::SeqConserved::Par seqPar_T; + seqPar_T.q = "Qpt_" + flavour[i] + "_5d"; + seqPar_T.action = "DWF_" + flavour[i]; + seqPar_T.tA = 0; + seqPar_T.tB = nt-1; + seqPar_T.curr_type = Current::Tadpole; + seqPar_T.mu_min = 0; + seqPar_T.mu_max = 3; + seqPar_T.mom = "0. 0. 0. 0."; + application.createModule("Qpt_" + flavour[i] + + "_seq_T", seqPar_T); + // seq propagator with tadpole insertion + MFermion::GaugeProp::Par quarkPar_seq_T; + quarkPar_seq_T.solver = "CG_" + flavour[i]; + quarkPar_seq_T.source = "Qpt_" + flavour[i] + "_seq_T"; + application.createModule("Qpt_" + flavour[i] + + "_seq_T" + flavour[i], + quarkPar_seq_T); + + + + //seq sources with conserved vector and photon insertion + MSource::SeqConserved::Par seqPar_V; + seqPar_V.q = "Qpt_" + flavour[i] + "_5d"; + seqPar_V.action = "DWF_" + flavour[i]; + seqPar_V.tA = 0; + seqPar_V.tB = nt-1; + seqPar_V.curr_type = Current::Vector; + seqPar_V.mu_min = 0; + seqPar_V.mu_max = 3; + seqPar_V.mom = "0. 0. 0. 0."; + seqPar_V.photon = "ph_field"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph", seqPar_V); + // seq propagator with conserved vector and photon insertion + MFermion::GaugeProp::Par quarkPar_seq_V; + quarkPar_seq_V.solver = "CG_" + flavour[i]; + quarkPar_seq_V.source = "Qpt_" + flavour[i] + "_seq_V_ph"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph_" + flavour[i], + quarkPar_seq_V); + + + + //double seq sources with conserved vector and photon insertion + //(for self energy) + MSource::SeqConserved::Par seqPar_VV; + seqPar_VV.q = "Qpt_" + flavour[i] + "_seq_V_ph_" + + flavour[i] + "_5d"; + seqPar_VV.action = "DWF_" + flavour[i]; + seqPar_VV.tA = 0; + seqPar_VV.tB = nt-1; + seqPar_VV.curr_type = Current::Vector; + seqPar_VV.mu_min = 0; + seqPar_VV.mu_max = 3; + seqPar_VV.mom = "0. 0. 0. 0."; + seqPar_VV.photon = "ph_field"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph" + flavour[i] + + "_seq_V_ph", seqPar_VV); + //double seq propagator with conserved vector and photon insertion + MFermion::GaugeProp::Par quarkPar_seq_VV; + quarkPar_seq_VV.solver = "CG_" + flavour[i]; + quarkPar_seq_VV.source = "Qpt_" + flavour[i] + "_seq_V_ph" + + flavour[i] + "_seq_V_ph"; + application.createModule("Qpt_" + flavour[i] + + "_seq_V_ph_" + flavour[i] + + "_seq_V_ph_" + flavour[i], + quarkPar_seq_VV); + + + + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + //2pt function contraction + MContraction::Meson::Par mesPar; + mesPar.output = "QED/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = ""; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + + + + //tadpole contraction + MContraction::Meson::Par mesPar_seq_T; + mesPar_seq_T.output = "QED/tadpole_pt_" + flavour[i] + "_T_" + + flavour[i] + "__" + flavour[j]; + mesPar_seq_T.q1 = "Qpt_" + flavour[i] + "_seq_T" + flavour[i]; + mesPar_seq_T.q2 = "Qpt_" + flavour[j]; + mesPar_seq_T.gammas = ""; + mesPar_seq_T.sink = "sink"; + application.createModule("meson_tadpole_pt_" + + flavour[i] + "_seq_T" + + flavour[i] + flavour[j], + mesPar_seq_T); + + + + //photon exchange contraction + MContraction::Meson::Par mesPar_seq_E; + mesPar_seq_E.output = "QED/exchange_pt_" + flavour[i] + "_V_ph_" + + flavour[i] + "__" + flavour[j] + "_V_ph_" + + flavour[j]; + mesPar_seq_E.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i]; + mesPar_seq_E.q2 = "Qpt_" + flavour[j] + "_seq_V_ph_" + flavour[j]; + mesPar_seq_E.gammas = ""; + mesPar_seq_E.sink = "sink"; + application.createModule("meson_exchange_pt_" + + flavour[i] + "_seq_V_ph_" + flavour[i] + + flavour[j] + "_seq_V_ph_" + flavour[j], + mesPar_seq_E); + + + + //self energy contraction + MContraction::Meson::Par mesPar_seq_S; + mesPar_seq_S.output = "QED/selfenergy_pt_" + flavour[i] + "_V_ph_" + + flavour[i] + "_V_ph_" + flavour[i] + "__" + + flavour[j]; + mesPar_seq_S.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i] + + "_seq_V_ph_" + flavour[i]; + mesPar_seq_S.q2 = "Qpt_" + flavour[j]; + mesPar_seq_S.gammas = ""; + mesPar_seq_S.sink = "sink"; + application.createModule("meson_selfenergy_pt_" + + flavour[i] + "_seq_V_ph_" + + flavour[i] + "_seq_V_ph_" + + flavour[i] + flavour[j], + mesPar_seq_S); + + } + + + + + // execution + application.saveParameterFile("QED.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} From d55212c998bba55b78bb676eb8a33f86fe6825cc Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Wed, 14 Feb 2018 10:45:18 +0000 Subject: [PATCH 7/8] restructure SeqConservedCurrent for DWF to need less memory --- lib/qcd/action/fermion/WilsonFermion5D.cc | 96 ++++++++++------------- 1 file changed, 43 insertions(+), 53 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 92280c13..6b830a03 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -883,8 +883,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); - PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), - tmp(FermionGrid()); + PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid()); Complex i(0.0, 1.0); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLs = q_in._grid->_rdimensions[0]; @@ -895,69 +894,60 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, 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*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) { - // 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); + bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); + bool tadpole_sign = (curr_type == Current::Tadpole); + bool switch_sgn = tadpole_sign || axial_sign; - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + + //forward direction: Need q(x + mu, s)*A(x) + ExtractSlice(tmp2, q_in, s, 0); //q(x,s) + tmp = Cshift(tmp2, mu, 1); //q(x+mu,s) + tmp2 = tmp*lattice_cmplx; //q(x+mu,s)*A(x) + + 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) { - 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; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteFwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, switch_sgn); } + } - // Repeat for backward direction. - t_mask = ((coords._odata[sU] >= (tmin + tshift)) && - (coords._odata[sU] <= (tmax + tshift))); + //backward direction: Need q(x - mu, s)*A(x-mu) + ExtractSlice(tmp2, q_in, s, 0); //q(x,s) + tmp = lattice_cmplx*tmp2; //q(x,s)*A(x) + tmp2 = Cshift(tmp, mu, -1); //q(x-mu,s)*A(x-mu,s) - //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 )); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + vInteger t_mask = ((coords._odata[sU] >= (tmin + tshift)) && + (coords._odata[sU] <= (tmax + tshift))); - timeSlices = Reduce(t_mask); + //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 )); - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) { - 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; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteBwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, axial_sign); } - } + } } } From 315a42843f654bf7e756bb30473ecfc129c11f26 Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Fri, 2 Mar 2018 11:47:38 +0000 Subject: [PATCH 8/8] changes requested for the pull request --- .../Hadrons/Modules/MSource/SeqConserved.hpp | 44 +++++---- extras/Hadrons/modules.inc | 2 +- lib/qcd/action/fermion/FermionOperator.h | 2 +- .../fermion/ImprovedStaggeredFermion.cc | 2 +- .../action/fermion/ImprovedStaggeredFermion.h | 2 +- .../fermion/ImprovedStaggeredFermion5D.cc | 2 +- .../fermion/ImprovedStaggeredFermion5D.h | 2 +- lib/qcd/action/fermion/WilsonFermion.cc | 71 +------------- lib/qcd/action/fermion/WilsonFermion.h | 4 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 94 +------------------ lib/qcd/action/fermion/WilsonFermion5D.h | 2 +- 11 files changed, 38 insertions(+), 189 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index 0fb23cc4..5c4a320b 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/extras/Hadrons/Modules/MSource/SeqConserved.hpp @@ -2,7 +2,7 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp +Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp Copyright (C) 2015-2018 @@ -35,8 +35,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include -#include - BEGIN_HADRONS_NAMESPACE /* @@ -100,6 +98,9 @@ protected: virtual void setup(void); // execution virtual void execute(void); +private: + bool SeqhasPhase_{false}; + std::string SeqmomphName_; }; MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); @@ -112,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); template TSeqConserved::TSeqConserved(const std::string name) : Module(name) +, SeqmomphName_ (name + "_Seqmomph") {} // dependencies/products /////////////////////////////////////////////////////// @@ -139,7 +141,7 @@ void TSeqConserved::setup(void) auto Ls_ = env().getObjectLs(par().action); envCreateLat(PropagatorField, getName(), Ls_); envTmpLat(PropagatorField, "src_tmp"); - envTmpLat(LatticeComplex, "mom_phase"); + envCacheLat(LatticeComplex, SeqmomphName_); envTmpLat(LatticeComplex, "coor"); envTmpLat(LatticeComplex, "latt_compl"); } @@ -152,37 +154,45 @@ void TSeqConserved::execute(void) { LOG(Message) << "Generating sequential source with conserved " << par().curr_type << " current at " - << "t = " << par().tA << std::endl; + << "t = " << par().tA << " summed over the indices " + << par().mu_min << " <= mu <= " << par().mu_max + << std::endl; } else { LOG(Message) << "Generating sequential source with conserved " << par().curr_type << " current for " << par().tA << " <= t <= " - << par().tB << std::endl; + << par().tB << " summed over the indices " + << par().mu_min << " <= mu <= " << par().mu_max + << 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); - 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; + auto &mom_phase = envGet(LatticeComplex, SeqmomphName_); + if (!SeqhasPhase_) + { + std::vector mom = strToVec(par().mom); + mom_phase = zero; + Complex i(0.0,1.0); + envGetTmp(LatticeComplex, coor); + 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); + SeqhasPhase_ = true; } - mom_phase = exp((Real)(2*M_PI)*i*mom_phase); - LOG(Message) << "Inserting momentum " << mom << std::endl; + LOG(Message) << "Inserting momentum " << strToVec(par().mom) << std::endl; diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 8bf22cc0..90602275 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -49,5 +49,5 @@ modules_hpp =\ Modules/MScalarSUN/TwoPoint.hpp \ Modules/MScalarSUN/TrPhi.hpp \ Modules/MIO/LoadNersc.hpp \ -Modules/MIO/LoadBinary.hpp + Modules/MIO/LoadBinary.hpp diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 94e065cc..9eb74baa 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -127,7 +127,7 @@ namespace Grid { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx)=0; + ComplexField &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 3a296e52..ed07fb5b 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -413,7 +413,7 @@ void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { assert(0); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 015eb6bb..d2426ef7 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -172,7 +172,7 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 3640a26b..3510b8f1 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -425,7 +425,7 @@ void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { assert(0); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 1c211938..e9cf2013 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -184,7 +184,7 @@ namespace QCD { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 5b14651b..2d9cf22d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -407,75 +407,6 @@ 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, @@ -484,7 +415,7 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 1ef6451e..ea25ed7f 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -174,13 +174,13 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { 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, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 8a795a1e..6f82aad2 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,98 +779,6 @@ 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, @@ -879,7 +787,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx) + ComplexField &lattice_cmplx) { conformable(q_in._grid, FermionGrid()); conformable(q_in._grid, q_out._grid); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index 6a7c6c7d..21da4c31 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -228,7 +228,7 @@ namespace QCD { unsigned int mu, unsigned int tmin, unsigned int tmax, - Lattice> &lattice_cmplx); + ComplexField &lattice_cmplx); }; }}